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CHAPTER  1 


INTRODUCTION 

1.1  General 

Analytical  solution  of  most  practical  nonlinear  structures  is  currently 
not  possible.  The  complex  nonlinear  partial  differential  equations  of 
equilibrium  arising  from  irregular  geometry,  large  displacements,  and  in- 
elastic material  behavior  are  intractable  with  current  analytical  techniques. 
In  many  instances,  the  governing  equations  cannot  even  be  explicitly  written 
in  closed  form. 

Consequently,  nonlinear  analysis  of  practical  structures  is  possible 
only  with  approximate  numerical  techniques.  Finite  differences  and  lumped 
parameter,  or  analog,  models  were  once  popular  methods  to  formulate  and 
solve  the  differential  equations.  Such  methods  are  well  suited  for 
two-dimensional  structures  of  regular  geometry.  However,  they  are  cumber- 
some for  irregular  geometries,  three  dimensional  problems,  and  transition 
regions  between  low  and  high  stress  gradients. 

The  finite  element  method  (FEM)  developed  during  the  past  two  decades 
provides  an  alternative  modeling  technique  for  structural  analysis.  A 
structure  is  first  idealized  by  an  assemblage  of  discrete  pieces,  or 
elements.  Within  each  element,  simple  expressions  approximate  the 
response  as  functions  of  unknown  quantities  specified  at  common  points 
(nodes)  between  elements.  Then,  instead  of  numerically  expressing  the 
partial  differential  equations  of  equilibrium,  variational  principles  are 
employed  to  generate  an  approximate  set  of  equilibrium  equations  directly 
in  terms  of  the  field  variables,  usually  nodal  displacements.  The  FEM 
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has  proven  highly  successful  in  linear  analysis  because  of  its  ability 
to  model  complex  geometries  and  variable  material  properties,  and  to 
combine  various  types  of  structural  elements.  The  FEM  is  currently 
proving  to  be  even  more  powerful  for  nonlinear  analysis  and  is  the  focal 
point  of  this  work. 

Development  of  the  FEM  coincided  with,  and  was  prompted  by,  concurrent 
advances  in  digital  computer  technology.  Without  modern  computers,  the 
FEM  would  be  relatively  useless  for  practical  analysis.  The  enormous 
number  of  computations  involved  in  setting  up  and  solving  the  governing 
equations  is  particularly  well  suited  for  computer  application.  However, 
finite  element  computer  software  has  received  surprisingly  little  formal 
attention  as  a research  topic.  Most  researchers  are  concerned  with 
theoretical  formulations  and  practical  applications  of  the  FEM.  Only 
recently  has  the  importance  of  "engineered"  finite  element  software  been 
realized  by  the  profession.  The  true  beauty  and  power  of  the  FEM  is  not 
realized  until  one  has  analyzed  an  extremely  complex  structure  with  a 
well  designed  user-oriented  finite  element  software  system. 

This  work  has  attempted  to  survey  the  state-of-the-art  in  nonlinear 
finite  element  technology  and  to  design  and  implement  a comprehensive, 
user-oriented  software  system  that  brings  these  techniques  to 
practicing  engineers,  researchers,  and  future  system  developers. 

1.2  Current  Nonlinear  Finite  Element  Systems 

Currently  available  finite  element  systems  can  be  separated  into 
three  distinct  categories  according  to  internal  structure,  ease  of  use, 
and  overall  generality.  Complete  reviews  of  these  and  many  special-purpose 


programs  for  linear  and  nonlinear  analysis  can  be  found  In  Ref.  [55]. 

A)  MARC  [41],  ANSYS  [68],  and  ADINA  — Both  MARC  and  ADINA  (a 
revised  NONSAP  [7])  resulted  from  university  research  projects  on  nonlinear 
finite  element  and  material  model  formulations.  aNSYS  was  developed  In 
private  Industry  for  linear  analysis  and  later  extended  to  include  a 
nonlinear  capability.  MARC  and  ANSYS  were  developed  during  the  late  1960's 
whereas  ADINA  Is  relatively  recent  (1975).  Free  versions  of  these  systems 
are  not  supported  and  generally  are  of  little  practical  value.  Proprietary 
versions  have  been  extensively  upgraded  and  are  now  In  widespread  use. 

Each  program  has  geometric  and  material  nonlinear  capability  with 
solid  mechanics  as  the  primary  application  area.  MARC  has  the  more 
extensive  nonlinear  material  modeling  capability.  Time  history  Integration 
and  modal  analysis  features  enable  nonlinear  dynamic  analysis.  Sulution 
procedures  are  based  on  Incremental  or  marching  techniques  that  trace 
the  approximate  load-displacement  path.  Recent  versions  of  these  programs 
support  the  general  Newton-Raphson  method  that  corrects  for  errors  intro- 
duced during  the  Incremental  solution. 

The  popularity  of  these  programs  Is  due  to  their  being  the  first  to 
offer  new  finite  elements,  nonlinear  material  models,  and  equation  solvers 
In  a usable  form.  However,  they  have  some  undesirable  characteristics. 

For  example.  Input  Is  generally  through  fixed-format  lines  of  data  In 
the  batch  mode  (the  most  recent  versions  are  beginning  to  support  a 
question-and-answer  Interactive  dialogue  type  of  input).  These  programs 
are  generally  not  complete.  Frequently,  users  must  develop  extensive 
pre  and  post  processing  packages  that  reduce  engineering  time  necessary 
to  specify  the  structural  model  and  Interpret  the  results. 
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The  Internal  program  and  data  structures  found  in  these  systems  are 
typical  of  software  generated  by  university  research  projects.  The 
simple  vector  and  matrix  data  structure  concepts  supported  directly  by 
FORTRAN  are  utilized.  No  integrated  data  management  or  virtual  memory 
techniques  are  employed.  The  advantages  of  this  approach  are  short 
development  time  and  low  development  cost. 

There  are  significant  disadvantages  associated  with  this  type  of 
finite  element  software.  From  the  user's  point  of  view,  these  programs 
appear  as  a hodgepodge  of  unrelated  packages.  Often  the  user  must 
execute  a series  of  programs  to  complete  a single  analysis.  Addition  of 
new  features  can  require  extensive  re-writes  of  the  system.  Extensions 
to  add  new  types  of  finite  elements  and  material  models  are  especially 
difficult  to  implement,  since  no  formal  communication  interfaces  exist 
between  program  processing  routines  and  element-material  model  routines. 
Similarly,  element  and  material  model  routines  are  not  logically  separate  — 
element  routines  invoke  material  model  routines  directly  thus  generating 
hidden  lower  level  interfaces.  Both  of  the  above  factors  greatly  compli- 
cate the  process  of  adding  new  elements  and  models.  These  disadvantages 
are  of  no  serious  consequence  for  the  casual  user,  but  researchers 
interested  in  augmenting  the  system  must  have  knowledge  of  the  internal 
organization,  existing  elements,  and  material  models  in  order  to  install 
new  elements  and  models.  In  addition,  developers  must  program  many  tasks 
such  as  input,  output,  and  memory  allocation,  that  are  not  directly 
related  to  the  element  or  material  model.  These  tasks  are  common  to  all 
elements  and  material  models  and  should  be  performed  by  the  system  itself. 
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B)  NASTRAN  [38],  ASKA  [5,66]  — These  systems  are  considerably 
more  complex  In  terms  of  Internal  organization  and  data  structure  than 
those  of  A)  and  were  developed  at  great  expense  (millions  of  dollars) 
through  government  funding.  They  were  originally  designed  for  static 
and  dynamic  analysis  of  very  large  linear  structures.  Neither  system 
has  a nonlinear  capability  comparable  to  the  research  oriented  programs 
of  A).  ASKA  Is  currently  undergoing  an  extensive  re-write  at  ISO;  at 
this  time  no  Information  Is  available  on  planned  nonlinear  capabilities. 
Proprietary  versions  of  the  NASTRAN  system  are  the  only  desirable  ones 
available  today.  These  have  updated  element  libraries  and  are  considerably 
more  efficient  than  the  public  domain  versions;  few  nonlinear  extensions 
have  been  Implemented.  Both  the  ASKA  and  NASTRAN  systems  still  receive 
considerable  financial  support  from  government  and  private  Industry 

The  popularity  of  NASTRAN  and  ASKA  stems  from  their  being  the  first 
to  solve  very  large  practical  structures.  They  are  still  not  user- 
oriented.  Input  Is  generally  fixed-format,  although  ASKA  now  does  have 
limited  free-form  input.  Numerous  pre  and  post  processing  packages  are 
available,  especially  for  the  NASTRAN  program. 

Developers  of  new  elements  and  material  models  face  the  same 
difficulties  as  with  the  programs  of  A).  System- element  Interfaces  are 
not  formalized  thereby  requiring  a knowledge  of  the  system  organization 
to  Install  new  elements  and  nonlinear  material  models. 

NASTRAN  and  ASKA  are  based  on  a concept  that  evolved  during  the 
mid  1960's  known  as  "prograitmlng  systems".  With  this  scheme,  a number 
of  processing  modules  are  designed  to  perform  specific  operations  on  data 
blocks  (NASTRAN)  or  hypermatrix  data  structures  (ASKA).  The  finite  element 
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system  then  consists  of  appropriate  sequences  of  calls  to  the  data  block 
and  hypermatrlx  processors.  NASTRAN  and  ASKA  have  pre-pro jrammed  sequences 
or  "formats"  for  standard  analysis  procedures;  I.e.,  static  linear  analysis, 
modal  analysis,  time  history  Integration,  etc. 

Processing  modules  request  data  through  a separate  data  rrianagement 
system.  NASTRAN  employs  DMAP-GINO,  a filing  system  that  operates  in  a 
simple  matrix  extraction  mode.  Data  management  functions  In  ASKA  are 
performed  with  DRS  (Data  Retrieval  System)  that  operates  on  hypermatrlx 
data  structures  stored  physically  as  orthogonal  lists. 

C)  STRUDL  [31],  FINITE  [34,35]  — These  two  systems  are  based  upon 
the  concept  of  Integrated  engineering  software  — a significantly  different 
approach  than  used  In  systems  of  B).  Both  use  sophisticated  supervisory 
systems  that  aid  engineers  In  the  development  of  complex,  user-oriented 
software  for  all  types  of  engineering  problems,  not  just  finite  element 
analysis.  Data  structures  considerably  more  general  than  hypermatrices 
are  supported  thereby  eliminating  many  of  the  difficulties  encountered  In 
systems  of  B). 

STRUDL  operates  under  the  ICES  supervisor  [61],  while  FINITE  operates 
under  the  POLO  II  supervisor  [32,33,37].  The  two  supervisors  perform 
essentially  the  same  tasks;  however,  the  philosophy  of  Implementation  Is 
radically  different. 

ICES  provides  software  development  support  through  ICETRAN,  a 
language  similar  to  FORTRAN  but  with  extended  data  structures  and  memory 
management.  ICETRAN  1s  converted  to  standard  FORTRAN  !5y  the  ICES 
precompiler  which  Inserts  numerous  statements  referring  to  ICES  run  time 


support  routines.  These  routines  are  merged  with  FORTRAN  emitted  by 
the  precompiler  to  form  a complete  program.  The  ICES  run  time  support 
routines  were  written  1n  machine  language  thereby  making  the  supervisor 
highly  machine  dependent. 

POLO,  developed  after  ICES  was  operational,  provides  more  compre- 
hensive data  management  capabilities  and  operates  In  a compiled  Interpre- 
tive mode.  The  run  time  support  routines  are  written  In  FORTRAN  thus 
making  POLO  considerably  more  machine  Independent  than  ICES.  Additional 
details  of  the  POLO  system  relative  to  nonlinear  finite  element  software 
are  given  In  Chapter  5. 

Emphasis  during  the  development  of  STRUDL  was  on  ease  of  defining 
the  structural  model  as  well  as  the  computational  aspects  of  analysis. 
However,  It  was  originally  designed  only  for  frame  analysis;  a finite 
element  capability  was  added  afterwards.  Some  versions  of  STRUDL  offer 
a limited  geometric  nonlinear  capability  for  buckling  analysis  of  frames. 
Nonlinear  material  behavior  Is  not  supported. 

The  FINITE  system  was  designed  to  have  the  same  ease  of  use  as 
STRUDL  but  with  greatly  extended  structural  modeling  capabilities  through 
multi-level  substructuring  and  static  condensation.  FINITE  Is  a true 
finite  element  system;  the  frame  analysis  features  so  popular  In  STRUDL 
are  Incorporated  In  the  same  manner  as  any  other  type  of  finite  element. 

A distinct  feature  of  the  FINITE  system  Involves  a formal  element 
''definition"  procedure.  The  system-element  connunl cation  Interface  Is 
standard  for  all  elements.  Knowledge  of  system  organization  Is  not 
required  to  Implement  new  elements.  In  addition,  FINITE  performs  many 
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features  cornnon  to  all  elements  that  must  be  programmed  by  developers 
in  the  systems  already  discussed.  As  part  of  this  study,  the  formal 

t » 

K f 

"definition"  procedure  was  extended  to  include  both  nonlinear  finite 
elements  and  material  models. 

1.3  Objectives  and  Scope 

The  objectives  of  the  research  reported  herein  are  threefold: 

1)  To  summarize  the  equations  of  nonlinear  continuum  mechanics 
in  matrix  form  suitable  for  application  to  finite  element 
analysis  and  to  determine  the  more  appropriate  formulation 
for  incorporating  large  displacement  effects; 

2)  To  design  a user-oriented,  general  static  nonlinear  finite 
element  system  that  includes  such  features  as  multi-level 
substructuring  and  condensation.  The  internal  organization, 
analytical  features,  and  user-system  interfaces  proposed  here 
should  serve  as  a model  for  software  to  satisfy  the  needs  of 
practicing  engineers  and  researchers.  Particular  emphasis  in 
this  work  is  placed  on  interaction  between  the  system  and 
finite  element  and  material  model  researchers. 

3)  To  implement  the  above  within  the  POLO-FINITE  structural  j 

mechanics  system  to  show  that  the  proposed  concepts  are  j 

feasible  and  efficient  for  solution  of  nonlinear  finite 
element  problems. 

Correspondingly,  this  report  is  divided  into  chapters  that  present 
the  techniques  used  to  achieve  the  objectives. 

Currently  available  numerical  formulations  for  nonlinear  finite 


element  analysis  are  reviewed  in  Chapters  2 and  3 to  establish  those 


procedures  applicable  to  a wide  class  of  problems.  Both  geometric  and 
material  non! Inearl  ties  are  considered.  Finite  strain  behavior  Is 
Incorporated  within  the  Lagranglan  formulation  adopted. 

Chapter  4 examines  the  various  aspects  of  finite  element  software 
that  directly  affect  users.  This  Includes  user-system  communication, 
structural  modeling  through  substructuring  and  condensation,  and  auto- 
matic analysis  restart. 

Chapter  5 Is  an  overview  of  the  nonlinear  POLO-FINITE  system  as 
designed  and  Implemented  during  this  work.  Particular  topics  addressed 
Include  subsystem  and  data  base  structure,  nonlinear  material  modeling 
for  the  analyst  and  material  behavior  researcher,  and  processors  for 
directing  the  nonlinear  analysis. 

The  solutions  to  several  example  problems  are  presented  In  Chapter 
6 to  Illustrate  the  applicability  of  FINITE's  nonlinear  capabilities. 

Chapter  7 briefly  sunmarlzes  this  research  effort. 
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CHAPTER  2 

NONLINEAR  STRUCTURAL  MECHANICS 


2,1  General 

All  rational  approaches  to  the  development  of  nonlinear  finite  element 
formulations  rely  upon  the  basic  principles  of  nonlinear  continuum  mechanics. 
This  chapter  presents  a concise  summary  of  the  salient  features  of  nonlinear 
continuum  mechanics  In  a form  suitable  for  direct  application  to  the  finite 
element  method. 

The  components  of  nonlinear  continuum  theory  necessary  to  develop  a 
finite  element  formulation  are: 

1)  Strain-displacement  relations; 

2)  Stress  definitions; 

3)  Constitutive  equations  relating  stress  to  strain; 

4)  Virtual  work  principles; 

Equations  describing  each  of  the  above  components  are  given  In  terms 
of  a general  three  dimensional  body  with  respect  to  rectangular  coordinates. 
Specialization  of  these  general  equations  result  In  explicit  formulations 
for  oriented  bodies  such  as  rods,  beams,  plates,  and  shells. 

No  attempt  has  been  made  here  to  derive  all  the  equations  from  first 
principles.  Rather  the  results  from  several  classical  treatises  are 
summarized  In  matrix  notation  familiar  to  most  structural  engineers.  For 
detailed  derivations,  the  reader  Is  referred  to  the  works  of  Murnaghan 
[45],  Novozhilov  [53],  and  Oden  [54]. 
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2.2  Coordinate  Systems 
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There  are  essentially  two  unique  ways  to  describe  the  deformation  of 
a body.  In  the  first,  all  quantities  of  Interest  are  expressed  as  func- 
tions of  the  rectangular  coordinates  of  the  body  In  Its  Initial  or  unde- 
formed state.  Such  a formulation  Is  termed  Lagranglan  or  Total  Lagranglan 
in  the  literature.  Strain  components  of  linear  elasticity,  given  In  terms 
of  displacement  derivatives  with  respect  to  the  Initial  coordinates,  are 
derived  from  a Lagranglan  approach. 

In  the  second  approach,  termed  Eulerlan,  convected  coordinates,  and 
Updated  Lagranglan,  rectangular  coordinates  are  continually  updated  to 
reflect  the  changing  configuration  of  the  body.  Subsequent  changes  of 
displacement,  strain,  and  stress  are  expressed  as  functions  of  the  Instan- 
taneous coordinates.  Linear  elasticity  makes  Implicit  use  of  the  Eulerlan 
system  to  derive  the  equilibrium  equations  by  considering  an  Infinitesimal 
volume  Isolated  from  the  deformed  body. 

If  only  infinitesimal  displacement  derivatives  are  present,  the  two 
approaches  are  essentially  identical.  However,  If  large  displacement 
gradients  occur,  resulting  In  excessive  rotations  and/or  deformation, 
governing  equations  derived  from  the  two  methods  are  different.  The 
choice  of  formulations  Is  a matter  of  convenience.  The  Lagranglan 
approach  Is  natural  for  expressing  the  deformed  geometry  and  strain  of 
a body  while  equilibrium  equations  are  simplest  in  the  Eulerlan  formulation. 
Geometric  transformations  relating  the  Initial  and  deformed  configurations 
of  a body  permit  the  use  of  both  formulations  In  a finite  element  analysis. 
These  transformations  are  derived  In  the  following  sections. 

Consider  an  arbitrary  three  dimensional  body  In  an  Initial  configura- 
tion (taken  as  undeformed  for  simplicity)  as  shown  In  Fig.  2.2.1  All 
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points  P(x,y,z)  are  located  relative  to  a rectangular  coordinate  system 
{X}  with  unit  vectors  {1}  directed  along  the  axes.  As  a result  of  an 
applied  loading,  each  point  P undergoes  a displacement  {U(x,y,z)}  with 
projections  onto  the  unit  vectors  {1}  denoted  by  [u,v,w].  The  new 
coordinates  of  points  In  the  body,  {X'}  are  given  by 


{X'}  * {X}  + {U} 


(2.2.1) 


where  the  terms  of  {X'}  are  functions  of  the  Initial  coordinates.  The 
Jacobian  matrix  of  the  functions  {X'}  relates  differentials  In  the  {X} 
and  {X ' } systems  by 


{dX'}  - [J]  {dX} 


(2.2.2) 


where. 


{dX}'  » [dx,  dy,  dz] 
{dX'}^  » [dx',  dy',  dz'] 


(2.2.3) 


(2.2.4) 


3(X'.Y'.Z' 

3(X.Y,Z) 


X'  X'  X' 
^x  > '^z 

Y'  Y'  Y' 
X y z 

Z'  Z'  V 
X y z 


(2.2.5) 


The  Jacobian  matrix  must  possess  a positive  determinant  for  all  points 
P(x,y,z).  A positive  determinant  assures  a unique  relationship  exists 
between  the  Initial  and  deformed  configuration,  I.e.,  two  points  In  the 
Initial  configuration  must  occupy  two  unique  positions  In  the  deformed 


configuration.  Furthermore,  the  [J]  matrix  defines  a linear  mapping  of 
lines,  surfaces,  and  volumes  In  the  Infinitesimal  neighborhood  of  each 
point  P.  Straight  lines  In  the  Initial  configuration  map  Into  straight 
lines  In  the  deformed  configuration.  Similarly,  planes  map  Into  planes 
and  parallelism  of  straight  lines  and  planes  Is  preserved. 

During  deformation,  coordinate  lines  originally  parallel  to  the  {X} 
axes  are  distorted  Into  a nonorthogonal  curvilinear  coordinate  system 
termed  {X}  as  shown  In  Fig.  2.2.2  for  a two  dimensional  body.  The  Jacobian 
matrix  contains  the  Information  describing  the  distortion.  Consider  a 
point  PQ(x,y)  and  a neighboring  point  P^(x  + dx,y  + dy)  that  enclose  an 
Infinitesimal  rectangle  of  area  dxdy.  After  displacements  have  occurred, 
the  rectangle  Is  distorted  Into  a parallelogram  with  new  corner  coordinates 
as  shown  In  the  figure.  In  particular,  sides  T dx  and  J dy  are  distorted 
Into  vectors  1 and  j.  From  simple  geometry,  the  projections  of  1 and  j 
onto  {1}  are  given  by  the  columns  of  [J].  The  area  of  the  parallelogram 
Is  given  by  1 x j or  simply  det[J]dxdy.  The  quantities  {dX'}  are  then  the 
projections  of  the  line  joining  the  points  Pq  and  P^  In  the  deformed 
configuration  on  unit  vectors  {1}.  Before  deformation  the  projections  were 


A differential  tetrahedron  In  the  Initial  configuration  with  centroid 
at  point  Pq  Is  shown  In  Fig.  2.2.3.  The  diagonal  face  has  area  dS  with 
an  outward  unit  normal  {n}.  Projections  of  dS  onto  the  coordinate  planes 


where  for  example  dA^^  Is  the  projection  of  dS  onto  the  YZ  plane.  After 
deformation  the  area  dS  Is  distorted  Into  dS'  with  a unit  outward  normal 
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{n'}.  The  projections  of  dS'  onto  the  coordinate  planes  are  denoted  {dA'}. 
By  considering  areas  initially  perpendicular  to  the  coordinate  axes  {X} 
and  following  the  same  procedure  as  above  for  the  two  dimensional  case, 
Murnaghan  [45]  showed  the  relationship  between  {dA'}  and  {dA}  can  be 
written  as 

{dA'}  = |J1  {dA}  (2.2.7) 

where  the  magnitude  of  area  change  is  related  to  both  |J|  and 
with  the  change  in  orientation  of  dS  provided  by  Similarly,  an 

area  in  the  deformed  configuration  with  projections  {dA'}  would  have 
projections  {dA}  given  by 

{dA}  = iJf^  [J]  {dA'}  (2.2.8) 

A differential  rectangular  volume,  dV,  in  the  initial  configura- 
tion becomes  a parallelepiped  with  volume  dV  in  the  deformed  configura- 
tion. The  ratio  of  deformed  and  undeformed  volumes  is  given  by 

dV  » |Jl  dV  (2.2.9) 

This  follows  from  a direct  application  of  the  triple  scalar  product  to 
compute  the  volume  of  a parallelepiped  using  the  columns  of  [J]  as  vectors 
tangent  to  each  edge. 

2.3  Strain-Displacement  Relations 

A major  source  of  nonlinear  behavior  arises  when  finite  rather  than 
infinitesimal,  displacement  gradients  and  rotations  occur  in  a body.  In 
this  section  the  relations  between  strain  and  displacement  for  arbitrarily 


large  displacement  gradients  are  presented  in  a form  suitable  for  finite 
element  analysis. 


2.3.1  Strain  Definitions 

The  definition  of  "engineering"  strain  was  first  proposed  by  Cauchy 
(1827)  and  is  of  fundamental  importance  in  linear  elasticity.  Cauchy 
defined  strain  in  terms  of  the  Initial  and  final  lengths  of  a fiber  in 
the  body  as 


where  e^.  denotes  Cauchy  strain.  and  are  the  initial  and  final 
lengths  of  a fiber. 

In  the  theory  of  finite  deformation,  two  definitions  of  direct 
strain  are  conmonly  employed.  The  first,  due  to  Green  (1839;,  expresses 
deformation  as  a function  of  the  initial  configuration  (Lagrangian  formula 
tion).  Green's  direct  strain  is  given  by 


The  factor  1/2  makes  Green  strain  equal  to  Cauchy  strain  for  small  length 
changes,  i.e.,  is  approximated  by  L^.  The  second  definition  of 
finite  direct  strain,  due  to  Almansi  (1911),  expresses  strain  as  a function 
of  the  deformed  configuration  (Eulerian  formulation).  Almansi 's  direct 
strain  is  given  by 


Again  the  factor  1/2  is  required  as  for  Green's  strain 
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2.3.2  Green  Strain  Components 

Components  of  direct  Green  strain  for  general  three  dimensional 
bodies  are  obtained  from  the  Jacobian  matrix  of  Eq.  2.2.5.  Consider  two 
points  a differential  distance  ds  apart  on  a space  curve  C before  deforma- 
tion as  shown  in  Fig.  2.3.1.  The  projections  of  ds  on  the  unit  vectors 
{i}  are  simply  {dX}.  The  magnitude  of  ds  is  the  positive  square  root  of 
{dX}^{dX}  . As  a result  of  deformation  the  curve  C is  distorted  into  a 
new  curve  C.  The  distance  between  the  two  points  after  deformation  is 
ds*  with  projections  {dX*}  on  unit  vectors  {i}.  Defining  ds*  as  the 
positive  square  root  of  {dX'}^{dX'},  the  following  relation  is  written 


(ds*)^  - (ds)^  = {dX*}''’  {dX*}  - {dX}’’’  {dX}  (2.3.4) 

2 2 

The  quantity  (ds*)  - (ds)  defines  the  deformation  in  the  infinitesimal 
neighborhood  of  the  point  P(x,y,z).  Writing  {dX*}  in  terms  of  {dX}  using 
[J]  yields 

(ds*)^  - (ds)^  = {dX}^  [J^J  - I]  {dX}  (2.3.5) 


If  the  distance  between  two  points  before  and  after  deformation  remains 
unaltered,  it  follows  that  [J^J]  = [I]  and  that  [J]  is  a simple  rotation 
matrix.  For  [J^J]  f [I],  the  difference  corresponds  to  that  part  of  the 
displacements  not  causing  rigid  body  motion  and  is  therefore  a measure  of 
the  deformation  at  a point.  In  view  of  the  above,  the  components  of  Green 
strain  are  thus  given  by 

[Eg]  = [J^O  - I]  1/2  (2.3.6) 

The  3x3  [sg]  matrix  contains  all  terms  necessary  to  describe  the  finite 
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deformation  of  a body  in  the  Laqrangian  formulation.  From  Eq.  2.3.6  the 
strain  components  are  seen  to  be  symmetric  and  are  usually  written  in 


vector  form  {cg}  for  finite  element  applications. 


Alternatively,  [cg]  can  be  written  in  terms  of  displacement  deriva- 


tives by  noting  that 


[J]  = [j]  + [I] 


(2.3.7) 


where 


[j]  = 


u u u 

X y z 


V V V 
X y z 


WWW 
X y z 


(2.3.8) 


Substituting  into  Eq.  2.3.6 


[Eg]  = [j  + + j^j]  1/2 


(2.3.9) 


Performing  the  operations  indicated  above,  the  components  of  Green  strain 


in  vector  form  {eg}  are 


- = u + i [u^  + v^  + w^] 
'X  X 2 '■  X X X-' 


e = V + [u^  + v^  + w^] 

y y 2 *■  y y y* 


= w + 1 [u?  + V^  + w^] 


^z  "z  2 ‘■“z  z "z 


(2.3.10) 


e =U  +V  + ruu  + VV  + wwl 
^xy  ”y  *x  ^'^x  y x y x y^ 


c =u  + w +ruu  + VV  + ww] 
XZ  Z X '■  Z X z X z x-^ 


'yz  " 'z  * «y  * ["y“2  * Vz  * Vz’ 


where  subscripts  imply  differentiation.  Clearly,  if  products  of  deriva- 
tives are  negligible,  the  components  of  Green  strain  simplify  to  those  of 
linear  elasticity. 

2.3.3  Almansi  Strain  Components 

Strains  for  the  Eulerian  formulation  defined  in  terms  of  the  instan- 
taneous coordinates  {X'}  are  those  due  to  Almansi.  The  inverse  of 
Eq.  2.2.2  provides  {dX}  in  terms  of  {dX'l 


{dX}  = [J]'^  {dX'l 


(2.3.11) 


Thus,  a line  segment  in  the  deformed  body  with  projections  {dX'}  has 
projections  {dX}  given  by  the  above  if  the  deformations  are  reversed. 
The  inverted  Jacobian  is  given  by 


m = [J] 


-1 


X,.  Xy.  X^. 

V- 

^x'  ^y'  ^z' 


or  in  terms  of  displacement  gradients. 


where. 


[j]  = [I]  - cn 


[H  = 


Ux.  Uy.  U^. 

h' 

Wx.  Wy.  W^. 


(2.3.12) 


(2.3.13) 


(2.3.14) 
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To  compute  [J],  the  displacements  (U)  are  added  to  the  coordinates  {X}, 
then  numerical  evaluation  of  [J]  and  [J]  are  Identical. 

Expressing  the  squared  change  In  length  of  a line  segment  as  a 
function  of  the  {X'}  coordinates.  Eq.  2.3.4  Is  rewritten  as 


(ds')^  - (ds)^  » {dX*}''^  [I  - f'S]  {dX’} 


from  which  the  components  of  Almansi  direct  strain  are 
[e,]  - [I  - fj]  1/2 
or  In  terms  of  displacement  gradients 


[e,]  • [J  + f - fH  1/2 


Almansi  strain  components  In  vector  form  {e_}  are 


Evl  “ “x*  ■ f 


Ey.  ” %•  " f 


e,,  - W^,  - J [u^,  + vj,  + wj,] 


S'y*  “ “y'  * ''x'  " ^•“x'^y'  ''x  y'  ”x'’^y'^ 


Ey.z'  “ * '^y*  “ t'*y''*z'  ''y'''z'  '*y''^2'^ 


(2.3.15) 


(2.3.16) 


(2.3.17) 


(2.3.18) 


If  the  products  of  differentials  are  negligible  and  the  displacements  are 
Infinitesimal,  the  above  strain  components  simplify  to  the  Cauchy  strains. 


2.3.4  Strain  Increments 

The  principle  of  virtual  work  for  large  deformations  derived  in  a 
subsequent  section  requires  an  expression  relating  strain  increments  to 
displacement  increments.  Increments  of  Green  strain  are  derived  below; 
Almansi  strain  increments  are  determined  with  an  identical  procedure. 

Consider  the  components  of  Green  strain  defined  in  Eq.  2.3.10.  For 
a variation  of  the  displacements  {6U},  the  variation  of  6t^  is  given  by 

^ I ^ ^ (2.3.19) 


which  can  also  be  written 


= (1  + Uj^)  (6Uj^)  + Vj^  (6Vj^)  + Wjj(6Wj^)  (2.3.20) 

The  coefficients  of  the  variational  terms  comprise  the  first  row  of  the 
Jacobian  matrix  of  Eq.  2.2.5.  For  infinitesimal  displacements  the 
derivatives  of  deformed  coordinates  {X'}  with  respect  to  undeformed 
coordinates  {X}  result  in  [J]  = [I]  with  the  above  variation  in  simpli- 
fying to  the  variation  of  Cauchy  strain.  Proceeding  as  above  for  each 
Green  strain  component,  the  resulting  variation  of  the  Green  strain  vector 
{6eg}  is  written  for  the  two  dimensional  case  as 


.}  = ' 6e..  > 


y , 

i'^xy' 


‘[B]  {6U} 


(2.3.21) 


where  [B]  is  a 3 x 2 differential  operator  defined  by 
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[B] 


¥»  ^ I Y*  — 
X 3x  X ax 


X'  ^ 


JL-  1 Y'  — 

'y  3y  y 3y 

yl  _L  I Y' 
y ax  ' y ax 


(2.3.22) 


v*  _5_  Y' 

X ay  I X ay  J 

Axlsynmetric  and  three  dimensional  expressions  for  the  [B]  matrix  are 
derived  In  a similar  manner. 


2.4  Stress 

In  general,  stress  provides  a measure  of  loading  Intensity  on  a 
specified  area  of  a body.  Relationships  among  stress  components  In  linear 
theory  are  derived  assuming  negligible  changes  in  geometry  due  to  displace 
ments.  In  such  cases,  stress  can  be  adequately  represented  as  functions 
of  the  Initial  coordinates  {X}.  However,  when  displacements  cause  large 
rotations  and/or  distortion  of  the  body,  several  definitions  of  stress 
arise  quite  naturally  as  functions  of  either  the  Initial  or  deformed 
coordinates.  Care  must  be  taken  to  assure  that  the  proper  stress  defini- 
tion Is  used  when  deriving  the  equilibrium  equations.  This  section 
examines  the  three  most  common  definitions  of  stress  proposed  In  finite 
deformation  theory,  namely,  Cauchy  or  Eulerlan  stress,  Lagranglan  or  1st 
Plola-KIrchoff  stress,  and  the  2nd  Plola-KIrchoff  stresses.  Equilibrium 
equations  for  the  body  are  derived  consistent  with  each  definition  of 


stress. 


2.4.1  Cauchy  Stress 

Consider  a differential  volume  of  a body  that  after  deformation  Is 
a simple  tetrahedron  as  shown  In  Fig.  2.4.1.  The  diagonal  face  has  area 
dS'  with  unit  outward  normal  {n'l.  Projections  of  the  vector  area  dS'{n'} 
on  the  coordinate  planes  are  {dA'l  as  shown  In  the  figure.  An  Infinitesimal 
force,  dT',  with  projections  {dT'}  along  the  unit  vectors  {1}  acts  over 
the  area  dS'.  This  force  arises  as  the  result  of  applied  loads  along  a 
bounding  surface  or  through  Interaction  with  adjacent  elements  of  the 
body.  The  3x3  matrix  of  stress  components  [o']  acting  on  the  projected 
areas  {dA'}  are  shown  In  the  figure.  Equilibrium  of  the  tetrahedron 


requires  that  stress  components  be  related  to  the  force  by 


{dT'}  = [oT  {n'}  dS' 


{dT'>  = [o']  (dA'l 


(2.4.1) 


(2.4.2) 


Elements  of  [o']  are  termed  Eulerlan,  Cauchy,  or  true  stresses  since 
they  reflect  actual  forces  acting  on  the  deformed  areas. 

Equilibrium  conditions  are  then  determined  by  equilibrating  the 
resultant  body  force  acting  on  any  arbitrary  volume,  v',  with  the  resultant 
of  tractions  applied  over  the  boundary.  S'  of  the  arbitrary  volume.  Integrat- 


ing Eq.  2.4.1  yields 


f (F'}  dV  = f [o'3{n': 
Jv  J<Z' 


(2.4.3) 


where  Integration  Is  with  respect  to  the  {X'}  coordinates  and  extends 
over  the  deformed  volume  and  surface  area.  The  vector  {F'l  contains 
components  of  applied  body  force  per  unit  deformed  volume.  Application 
of  the  divergence  theorem  to  the  right  side  of  Eq.  2.4.3  converts  the 


23 

surface  Integral  to  a volume  Integral  with  the  resulting  equilibrium 
conditions  given  as 

► 

[{F'}  + dIVy,  [aQ]  dV  » {0}  (2.4.4) 

JV  j 

The  above  can  hold  for  any  arbitrary  portion  of  the  body  only  If  the  Integrand 
vanishes  everywhere  In  the  body,  resulting  In  the  equilibrium  equations  In 
terms  of  stress  components 

d1vj(,  [o']  + {F'}  = {0}  (2.4.5) 

Similarly,  equating  the  total  moments  of  body  forces  and  surface  forces 
about  the  origin  shows  that  = oj^.  At  first  these  equilibrium  equa- 

tions appear  Identical  to  those  derived  In  linear  elasticity.  The  differ- 
ence Is  that  the  Eulerlan  stress  components  are  functions  of  the  coordinates 
{X'}  not  {X)  and  that  differentiation  is  with  respect  to  {X'}.  However, 
derivatives  with  respect  to  {X'}  cannot  be  obtained  until  displacements 
are  known.  This  difficulty  is  overcome  by  transforming  the  equilibrium 
equations  Into  functions  of  the  {X}  coordinates  as  shown  In  the  following 
sections. 

2.4.2  Lagrange  and  Plola-KIrchoff  Stresses 

The  limits  on  Integrals  of  Eq.  2.4.3  expressing  equilibrium  conditions 
extend  over  the  deformed  volume  and  surface  area  of  the  body.  Limits  on 
the  Integrals  are  transformed  by  v/rlting 

{F'XdV’/dV)  dV  - [o']{n'}(dS7dS)  dS  (2.4.6) 

Jv  Js 

such  that  the  Integrals  are  now  evaluated  over  the  Initial  configuration 

i 

of  the  body.  Substituting  for  the  components  of  {dA'}  In  terms  of  (dA)  i 

from  Eq.  2.2.7  and  for  dV  from  Eq.  2.2.9,  the  above  equation  can  be  i 


wri tten 
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{F'}1J|  dV  = [o']  1J|  [j''’]‘^{n}  dS  (2.4.7) 

V Js 

where  the  terms  of  the  3x3  matrix  defined  as 

[0^]  = [o']  |J|  [/]•■'  (2.4.8) 

are  called  the  Lagrange  or  1st  Piola-Kirchoff  stress  components.  Lagrange 
stress  components  act  in  directions  corresponding  to  unit  vectors  {i} 
over  the  initial  area  dS  and  equilibrate  the  same  differential  force  {dT'} 
as  do  Eulerian  stresses  acting  on  the  deformed  area  dS'.  The  same  result 
is  obtained  by  first  transforming  from  dS'  to  dS  in  Eq.  2.4.1  then 
expressing  the  equilibrium  condition. 

The  Jacobian  matrix  [J]  is  in  general  not  symmetric  thus  rendering 
the  components  of  Lagrange  stress  nonsymmetric  as  given  by  the  transforma- 
tion in  Eq.  2.4.8.  This  proves  to  be  inconvenient  as  both  Green  and 
Almansi  strain  components  are  symmetric  thus  necessitating  a nonsymmetric 
constitutive  relation  even  for  isotropic  materials. 

An  alternative  definition  of  stress  is  derived  by  noting  that  Lagrange 
stresses  acting  on  undeformed  areas  equilibrate  the  same  force  {dT'} 
acting  on  the  deformed  area.  However,  a force,  (dT),  acting  on  the 
undeformed  area  can  be  defined  by 

{dT}  = [J]'^{dT'}  (2.4.9) 

This  expression  states  that  a force  vector  with  components  {dT}  acting 
on  the  body  before  deformation  has  components  {dT'}  after  deformation 
which  is  analogous  to  the  changing  direction  of  a vector  connecting  two 
points  before  and  after  deformation.  Substituting  for  {dT'}  into  Eq.  2.4.1 
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3X  * 3y  j2*  * 

2.5  Principle  of  Virtual  Work 

The  principle  of  virtual  work,  or  more  precisely  the  principle  of 
virtual  displacements,  provides  an  alternative  but  equivalent  statement 
of  the  equilibrium  conditions  for  a body.  In  the  finite  element  method 
this  principle  is  used  to  derive  a finite  number  of  nonlinear  algebraic 
equations  that  approximate  the  nonlinear  partial  differential  equations 
of  equilibrium  in  Eqs.  2.4.5  and  2.4.14. 

Virtual  work  principles  are  customarily  derived  first  for  a particle, 
a collection  of  particles,  and  finally  for  a continuum.  Highlights  of  the 
derivation  for  the  Lagrangian  system  are  presented  here  since  it  is  more 
complex  than  the  Eulerian  derivation.  The  algebraic  operations  are  quite 
lengthy  but  relatively  straightforward.  Equilibrium  equations  of  2.45  are 
multiplied  by  arbitrary  virtual  displacements,{6U},  and  integrated  over  the 
deformed  volume.  To  this  is  added  the  integral  over  the  deformed  surface 
the  product  of  virtual  displacements  times  {T'}-{T'},  where  {T'}  are  components 
of  specified  surface  traction  and  {T'}  = [a']{n'}.  The  virtual  displacements 
are  considered  functions  of  the  initial  coordinates  {X}  and  satisfy  all 
geometric  boundary  conditions  imposed  upon  the  body.  The  integration  limits 
are  transformed  to  the  undeformed  configuration  by  introducing  the  Lagrange 
stress  components  as  in  Eq.  2.47.  Green's  theorem  is  then  applied  to  the 
surface  integral  transforming  it  into  an  equivalent  integral  over 
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the  undeformed  volume.  All  Lagrange  stress  components  multiplied 
by  variations  of  displacement  gradients  are  grouped  together, 
then  the  Lagrange  stress  components  are  expressed  in  terms  of 
2nd  Piola-Kirchoff  stresses  defined  by  Eq.  2.4.11.  The  coefficients 
multiplying  each  component  of  2nd  Piola-Kirchoff  stress  are  seen  to  be 
precisely  the  variations  of  Green  strain  components  defined  in  Eq.  2.3.22. 
The  remaining  terms  containing  derivatives  of  Lagrange  stresses  multiplied 
by  virtual  displacements  are  re-expressed  as  a surface  integral  giving 
the  virtual  work  of  all  surface  tractions  applied  over  the  body.  The 
resulting  integrals  provide  the  equilibrium  equations 


{6U}^{F}  dv  = {6e^}^  {o}  dV - 

I V 9 


{6U}'  {T}  dS 


(2.5.1) 


The  left  most  integral  is  the  virtual  work  of  all  body  forces  acting 
through  the  virtual  displacements.  The  middle  integral  represents  the 
internal  virtual  work  and  the  rightmost  integral  is  the  work  of  applied 
surface  tractions,  {T},  specified  in  directions  of  unit  vectors  {1}  with 
intensities  always  based  on  the  initial  area.  Defining  the  work  of  body 
forces  as  the  internal  work  6W^  as  the  negative  of  the  second 

integral, and  the  work  of  surface  tractions  as  6W^,  the  equilibrium  equations 


(2.5.2) 


As  will  be  shown  in  Chapter  3,  Eq.  2.5.1  is  directly  applicable  in  the 
nonlinear  finite  element  method. 
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The  virtual  work  equations  for  the  Eulerian  formulation  are  identical 
to  those  above  with  the  appropriate  changes  in  components 


{6U}^{F'}  dV  = {o'}  dV  - 

y I y ' ® 


{6U}^  {T'}  dS' 


(2.5.3) 


Variations  of  Almansi  strain  replace  Green  strain,  Cauchy  stresses  replace 
2nd  Piola-Kirchoff  stresses,  surface  tractions  are  for  deformed  areas  and 
body  forces  are  per  unit  deformed  volume. 

2.6  Stress-Strain  Relations 

The  finite  element  method  has  provided  new  impetus  for  research 
efforts  into  the  understanding  of  nonlinear  material  behavior.  Before 
finite  element  technology  was  available,  nonlinear  analysis  was  complicated 
by  material  behavior  as  well  as  large  deformations  and  irregular  boundary 
conditions.  The  finite  element  method  can  approximate  the  latter  two 
factors  to  within  any  accuracy  desired  but  must  still  rely  upon  theoretical 
models  of  material  behavior  for  stress-strain  relations. 

The  purpose  of  a material  model  is  to  predict  the  stress  at  a point 
in  the  body  given  the  strain.  This  may  be  accomplished  in  either  total 


or  incremental  form  as 


{o}  = [Dj]  {£> 


(2.6.1) 


{o}  = E{Aa}  = E[D^]{Ae}  (2.6.2) 

where  [D^]  is  termed  the  secant  modulus  matrix  and  [D^]  the  tangent  modulus. 


Strains  are  either  Green  or  Almansi  with  corresponding  2nd  Piola-Kirchoff 
stresses  or  Cauchy  stresses.  For  linear  elastic,  small  displacement 
analysis  [D^]  = [0^]  and  each  matrix  contains  the  elastic  constants 
referred  to  the  initial  configuration. 

When  stresses  exceed  the  proportional  limit  or  strains  become  finite 
in  magnitude,  [D^]  and  [D^]  are  no  longer  ^iven  by  the  elastic  constants 
and  a nonlinear  [D]  matrix  dependent  upon  current  strain,  previous  loading 
history,  etc.  is  required.  Three  basic  situations  arise: 

1)  Infinitesimal  strain  magnitudes  with  stresses  that  exceed 
the  proportional  limit  of  the  material.  Behavior  is  often 
adequately  predicted  by  classical  plasticity  and  fracture 
(cracking)  theories. 

2)  Strains  have  finite  magnitude  but  the  material  is  highly 
elastic,  for  example  rubber.  Mooney-Rivilin  models  have  been 
proposed  with  terms  of  [D^]  given  by  differentiation  of  the 
strain  energy  density  function  (determined  experimentally)  with 
respect  to  Green  strain  components. 

3)  Finite  strains  in  materials  that  exhibit  yielding  behavior  thus 
requiring  the  application  of  large  strain  plasticity  theory. 

Fortunately,  most  nonlinear  problems  of  interest  to  structural 
engineers  involve  the  first  type  of  behavior  listed  above  with  materials 
such  as  concrete  and  steel.  Structures  having  this  characteristic  include 
frames,  plates,  shells,  and  most  solid  bodies.  Frames,  plates  and  shells 
may  have  large  rotations  under  loading  but  small  strains. 
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A general  approach  to  nonlinear  finite  element  analysis  necessitates 
an  incremental  technique  for  solving  the  equilibrium  equations;  thus,  only 
incremental  forms  of  stress-strain  relations  of  Eq.  2.6.2  are  required. 

In  the  Lagrangian  formulation,  increments  of  2nd  Piola-Kirchoff  stress 
can  be  related  directly  to  increments  of  Green  strain  through  the  tangent 
modulus  matrix  for  the  material.  No  coordinate  transformations  are 
necessary  as  stresses  are  always  based  on  undeformed  areas  and  the  direc- 
tion of  strain  and  stress  components  coincide  with  unit  vectors  {i}. 

However,  Cauchy  or  Eulerian  stresses  based  on  projections  of  deformed 
areas  and  Almansi  strain  based  on  deformed  coordinates  cannot  be  directly 
related  as  the  deformed  area  and  coordinates  change  during  the  increment. 

For  an  increment  consisting  of  pure  rigid  body  rotation,  projections  of 
the  deformed  area  onto  coordinate  planes  change  thus  changing  the  Cauchy 
stresses.  A technique  to  account  for  rotations  during  an  increment  was 
introduced  by  Jaumann  who  defined  a stress  increment  {AojJ=[Dj]{Aeg} 
with  the  increment  of  Eulerian  stress  then  given  by 

{Aa  } = (Aoj)  + [aTJ  (2.6.3) 

where  [aT^]  is  the  matrix  of  incremental  rigid  body  rotations  due  to  the 
displacement  increment.  Thus,  the  change  of  stress  is  determined  by  a 
combination  of  rigid  body  rotation  and  additional  straining.  The  [Dj] 
matrix  contains  the  usual  elastic  constants  or  terms  resulting  from  classical 
plasticity.  The  terms  of  [aT^]  are  given  in  Ref.  [49].  Eq.  2.6.3  also 
shows  that  before  Eulerian  stress  increments  can  be  combined  they  must 
first  be  rotated  to  a common  coordinate  system  then  summed  and  transformed 
back  to  final  deformed  coordinates.  This  is  usually  accomplished  by 
transforming  to  2nd  Piola-Kirchoff  stress  increments  via  Eq.  2.4.11. 
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An  additional  complication  with  the  Eulerian  formulation  involves 
materials  that  are  initially  anisotropic.  As  the  body  deforms,  principal 
material  property  directions  constantly  change:  thus,  the  [D]  matrix  must 
be  transformed  at  each  increment  in  the  analysis.  This  problem  does  not 
occur  in  the  Lagrangian  formulation. 

2.7  Summary 

In  this  chapter  the  basic  theory  of  nonlinear  continuum  mechanics 
has  been  presented  in  a form  suitable  for  finite  element  application.  As 
shown,  two  uniquely  different  approaches,  Lagrangian  and  Eulerian,  are 
available  to  formulate  nonlinear  structural  analyses. 

Early  attempts  at  nonlinear  finite  element  analysis  employed  the 
Eulerian  formulation  as  a natural  extension  of  existing  linear  analysis 
systems.  Coordinates  of  node  points  were  simply  updated  by  the  incremental 
displacements  and  another  linear  analysis  performed.  However,  as  seen  in 
this  chapter,  such  an  approach  requires  continual  transformations  of 
coordinates,  strains,  stresses,  and  material  properties.  This  factor 
severely  limits  the  efficiency  of  the  Eulerian  approach  for  large  scale 
structural  analyses.  It  has  traditionally  been  most  successful  for 
applications  in  which  bodies  suffer  extreme  deformations  under  loading, 
such  as  a viscous  fluid. 

The  Lagrangian  formulation  has  none  of  these  disadvantages  and 
considerably  simplifies  the  overall  analysis.  It  is  applicable  for  all 
common  structures  that  do  not  exhibit  severe  distortions  under  loading. 
Therefore,  the  Lagrangian  formulation  was  adopted  in  this  study  for  imple- 
mentation of  a general  nonlinear  finite  element  system  for  structural  analysis. 
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CHAPTER  3 


NONLINEAR  ASPECTS  OF  THE  FINITE  ELEMENT  fCTHOD 


3.1  General 

The  finite  element  method  (FEM)  has  gained  widespread  acceptance 
among  engineers  for  the  analysis  of  linear  structures.  The  development 
of  elements  capable  of  modeling  most  structural  components  combined  with 
the  intuitive  and  mathematical  bases  of  the  method  has  contributed  to  its 
popularity.  Finite  element  techniques  have  also  proven  invaluable  in 
nonlinear  applications.  Nonlinear  analysis  of  practical  structures  is 
achieved  by  approximating  the  complex  nonlinear  partial  differential 
equations  of  continuum  mechanics  with  a set  of  nonlinear  algebraic  equa- 
tions. 

The  following  three  classes  of  nonlinear  structural  response  are 
considered  in  this  work: 

1)  Small  displacements,  small  strains  with  linear  and 
nonlinear  material  properties; 

2)  Large  displacements  resulting  in  significant  rotations  but 
small  strains  with  linear  and  nonlinear  material  properties. 

3)  Large  displacements  with  associated  finite  strain  magnitudes 
also  with  linear  and  nonlinear  material  properties. 

The  first  category  has  received  the  most  attention  of  finite  element 
researchers  and  includes  ordinary  linear  structures  and  solid  mechanics 
problems  with  material  nonlinearities.  Elastic  and  inelastic  stability 
analyses  of  frames,  plates  and  shells  fall  into  the  second  category. 
Problems  in  which  finite  strains  occur  are  generally  the  concern  of 
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engineers  working  In  non-structural  applications  and  Involve  materials 
such  as  rubber  and  plastic.  Emphasis  In  this  study  was  placed  on  the 
first  and  second  classes  of  structures  for  which  some  constitutive  relations 
are  available  for  common  construction  materials.  However,  the  computer 
system  discussed  In  Chapter  5 based  on  the  formulation  presented  here.  Is 
capable  of  solving  finite  strain  problems  If  appropriate  constitutive 
relations  are  provided. 

This  chapter  briefly  examines  the  finite  element  method  for  nonlinear 
analysis  within  the  Lagranglan  formulation.  Full  details  of -the  method 
for  the  first  class  of  problems  listed  above  are  available  In  most  standard 
texts  D3, 15,73].  The  discussion  here  concentrates  on  extensions  of  the 
basic  formulation  necessary  to  Incorporate  large  displacement  effects. 
Details  of  these  techniques  are  not  generally  available  In  the  literature. 
The  response  of  Isolated  Individual  elements  Is  discussed  prior  to  pro- 
cedures for  generating  and  solving  the  nonlinear  equilibrium  equations. 
Special  consideration  Is  given  to  the  residual  load  concept  and  the  appli- 
cation of  substructuring  with  static  condensation. 

3.2  Individual  Elements 

Once  a structure  is  partitioned  Into  a mesh  of  finite  elements. 
Individual  elements  may  be  Isolated  and  studied  Independently.  Nonlinear 
equations  describing  the  response  of  a single  element  are  considered  In 
this  section.  Discussion  Is  limited  to  formulations  based  on  assumed 
displacement  fields  as  these  elements  are  most  commonly  employed  In 
practical  analysis.  Strain-displacement  relations  and  element  equilibrium 
equations  are  derived  following  a discussion  of  Interpolating  functions. 
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A general  three  dimensional  element  in  rectangular  coordinates,  as  shown  in 
Fig.  3.2.1,  is  employed  to  illustrate  specific  terms  of  the  matrices  in- 
volved. 


3.2.1  Interpolation  Functions 

The  conversion  of  discrete  displacement  component  values  specified 
at  node  points  into  continuous  analytic  functions  of  the  element  local 
coordinates  constitutes  the  basic  step  of  the  displacement  formulation. 
Displacement  components  (degrees  of  freedom,  dof)  at  each  node  are  listed 
in  an  M X 1 vector,  denoted  (1)^},  where  M is  the  total  number  of  dof 
for  the  element  nodes.  Displacements  at  points  P(x,y,2)  within  the  element 
as  functions  of  local  element  coordinates  are  given  by 

nn 

{U(x,y,z)}  = E [N^(x,y,z)]  {Ug}^.  (3.2.1) 

where  nn  is  the  number  of  element  nodes,  [N^]  is  a matrix  of  interpolating 
or  "shape"  functions  for  the  i^^  node,  and  {Ug}^  are  displacement  components 
for  the  i^^  node.  For  a three  dimensional  element  with  dof  U,  V,  and  W 
at  each  node,  the  above  equation  is 


1 U(x,y,z) 

nn 

( V(x,y,z) 

\ = r 

( , 

W(x,y,z) 

1 

. 

N.  0 0 


0 

0 0 


N.  0 

Ni 


(3.2.2) 


Each  component  of  displacement  above  is  interpolated  over  the  element 
with  identical  shape  functions.  In  general,  different  dof  can  be  defined 
at  each  node  thereby  forming  a variable  dof  element.  Such  an  element  proves 
useful  to  transition  between  two  different  classes  of  elements,  for  example, 
shells  (6  dof/node)  and  three  dimensional  solids  (3  dof/node). 
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• Other  quantities  for  the  element  may  also  be  Interpolated  with  shape 

f. 

functions.  For  example,  the  thickness  of  planar  elements,  varying  material 
properties,  temperature,  etc.  are  easily  converted  to  functions  of  element 
coordinates  by  Interpolation  from  specified  nodal  values. 

3.2.2  Strain-Displacement  Relations 

Before  the  element  equilibrium  equations  can  be  derived.  It  Is  first 
necessary  to  determine  the  variation  of  fireen  strain  corresponding  to  a 

I 

I 

i variation  of  nodal  displacements.  From  Eq.  2.3.21  the  required  relation 

I 

can  be  written 

I 

^ {fiEg}  = [B]  {6U}  (3.2.3) 

j where  {fiCg}  Is  a 6 x 1 vector  containing  the  variation  of  Green  strain 

components  as  functions  of  the  undeformed  element  coordinates.  [B]  Is 
I the  6x3  differential  operator  matrix  defined  In  Eq.  2.3.22.  {6U}  Is 


the  3 X 1 vector  of  local  element  displacement  variations  defined  In  terms 
of  nodal  displacement  variations  by 

nn 

{6U}  = E [Nj  {6U^}.  (3.2.4) 

1 1 e 1 

Derivatives  of  the  deformed  coordinate  system  {X'}  appearing  in  terms  of 
[B]  are  obtained  by  differentiating 

nn 

{X'}  = {X}  + E [N.]  {U^}^  (3.2.5) 

1 ^ ® ^ 

Similarly,  derivatives  of  displacements  are  given  by  differentiating 


Eq.  3.2.1  Substituting  derivatives  In  terms  of  shape  functions  Into 
Eq.  2.3.22  yields  the  variation  of  Green  strain  In  terms  of  finite  element 
quantities  for  the  three  dimensional  case  as 
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nn 

{«£g}  = E [B^]  {6Ug}. 

where  the  terms  of  [B^]  are 
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(3.2.6) 


(3.2.7) 


Derivatives  of  shape  functions  in  the  equation  refer  to  the  function  for 
the  i^^  node,  N..  For  small  displacement,  linear  analysis,  X'  = Y'  = Z'  = 1; 

I A jr 

other  derivatives  of  {X'}  are  zero.  [B^.]  then  simplifies  to  the  well 
known  matrix  employed  in  linear  analysis.  However,  in  large  displacement 
analysis,  the  [B]  matrix  is  a function  of  the  nodal  displacements  because 
derivatives  of  deformed  coordinates  do  not  vanish.  By  imposing  the 
Kirchoff  plane  sections  hypothesis,  specific  forms  of  the  [B]  matrix  are 
obtained  for  frames,  plates,  and  shells. 
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An  expression  for  the  total  Green  strain  follows  readily  from  terms 
of  the  [B]  matrix.  First,  separate  [B^]  Into  the  sum  of  two  matrices, 
[B^]  and  [B^*"],  such  that  [B^]  represents  the  usual  linear  matrix  with 
[B"  ] containing  the  nonlinear  terms.  A simple  re-examination  of  the 
expanded  total  Green  strain  components,  Eq.  2.3.10,  shows  that  the  total 
strain  Is  given  by 


(tg)  " £ [[B^]  + i [b'}'-]]  (U^), 


(3.2.8) 


3.2.3  Element  Loads 


Loads  applied  over  an  element  are  expressed  as  functions  of  local 
element  coordinates  by  Interpolation  of  nodal  values.  Intensities  of 
body  force  components,  {F},  are  given  by 


{F}  = E [NJ  {F^}. 
1 1 el 


where  nodal  values  are  force  per  unit  undeformed  volume. 


(3.2.9) 


Components  of  a traction  applied  over  the  element  boundaries  are 


Interpolated  from  traction  Intensities  at  nodes  such  that 


fill 

{T}  = E [N^]  {Tg}^ 


(3.2.10) 


where  Is  the  3 x 1 vector  of  traction  components  at  the  node 

referred  to  Initially  undeformed  areas  and  element  local  coordinate 
directions. 


3.2.4  Element  Virtual  Work  Equations 

All  components  required  to  write  the  virtual  work  equation  In  the 
Lagranglan  formulation,  Eq.  2.5.1,  are  now  available  In  terms  of  finite 


I 


element  quantities.  By  substituting  Eqs.  3.2.4  and  3.2.9  into  Eq.  2.5.1, 
the  virtual  work  of  body  forces  is  seen  to  be 


[N^]  {Fg}^.  dV  (3.2.11) 

A simpler  form  of  this  equation  is  obtained  by  using  an  implied 
summation  over  the  element  nodes  and  by  removing  nodal  quantities  inde- 
pendent of  the  coordinates  under  the  integral.  Thus,  Eq.  3.2.11  simpli- 
fies to 
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bf 


nn 
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= {6U}'  [ 


[N]'  [N]  dV]  {F} 


(3.2.12) 


where  integration  is  over  the  initial  volume  of  the  element. 

Similarly,  the  virtual  work  of  applied  element  surface  tractions 
is  given  by 


[N]'  [N]  dS]  {T} 


(3.2.13) 


with  integration  performed  over  the  undeformed  element  surface. 
The  total  virtual  work  of  applied  element  loads  is  then 


= {6U}-  {P} 


where. 


{P}  = [ 


[N]'  [N]  dV]  {F}  + [ 


(3.2.14) 


[N]^  [N]  dS]  {T}  (3.2.1£() 


{PI  contains  the  equivalent  (consistent)  nodal  loads  for  the  element. 
Nonlinear! ties  enter  Eq.  3.2.15  in  large  displacement  analysis 


when  significant  differences  exist  between  the  initial  and  deformed 
element  geometry.  In  such  cases,  body  force  and  surface  traction  components 


i . ; 


at  nodes  specified  In  the  deformed  geometry  must  be  transformed  to  the 
Initial  geometry  before  evaluation  of  Eq.  3.2.15. 


The  Internal  virtual  work  Is  simply 


il  [Bf 

Jv 


(o)  dV) 


(3.2.16) 


In  finite  element  terms.  This  can  also  be  written  as 


6Wint  = -{6U>  {IF} 


(3.2.17) 


where  (IF)  Is  often  called  the  element  Internal  nodal  force  vector. 

Terms  of  {IF}  are  components  of  force.  In  the  generalized  sense,  exerted 
on  the  element  by  the  nodes  due  to  the  element's  state  of  deformation.  Non- 
linearities  enter  through  the  [B]  matrix  for  large  displacements  and  through 
{a}  for  nonlinear  material  behavior. 


3.3  Structure  Equilibrium  Equations 

The  structural  equilibrium  equations  In  terms  of  nodal  forces  result 
from  a straightforward  application  of  statics  at  each  node.  In  summing 
forces  at  each  structural  node,  element  Internal  forces  and  equivalent 
nodal  loads  require  rotation  from  element  local  to  structural  global 
coordinate  frames  and  proper  placement  In  the  structure  equilibrium 
equations.  This  Is  accomplished  In  the  usual  way  with  an  element  rotation 
matrix,  [x^],  and  the  Boolean  connectivity  matrix,  [L^],  such  that  the 
structural  equilibrium  equations  are  generated  symbolically  as 

(R}  = I [L^f  [A^flp}^  - t [L^f  [X^f  {IF}^  (3.3.1) 


or  more  simply. 


{R}  =•  {Pj}  - {IFj} 


(3.3.2) 
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where  {P^}  and  {IF^}  are  the  total  structure  applied  load  vector  and 
total  internal  force  vector  respectively.  The  summation  over  all  elements 
also  includes  substructures  in  addition  to  simple  finite  elements  with- 
out any  changes  in  notation.  The  difference  between  applied  load  and 
internal  force  vectors,  termed  the  residual  nodal  loads  {R},  must  vanish 
for  equilibrium. 

Eq.  3.3.1  shows  a major  computational  advantage  of  the  Lagrangian 
formulation.  The  rotation  matrix,  [x^],  for  each  element  is  independent 
of  the  nodal  displacements.  In  the  Eulerian  approach,  every  modification 
of  nodal  displacements  necessitates  recomputation  of  the  rotation  matrices. 

3.4  Solution  of  Equilibrium  Equations 

Each  term  of  the  residual  nodal  load  vector  {R}  is  indirectly  a 
nonlinear  algebraic  function  of  the  ncJal  displacement  components  {U}. 

A set  of  nodal  displacements  corresponding  to  an  equilibrium  configura- 


tion satisfies  the  relation 


{R({U})}  = {Oi 


(3.4.1) 


Solution  of  the  nonlinear  finite  element  problem  requires  solution  of 
'n'  simultaneous  nonlinear  algebraic  equations  where  'n'  is  the  total 
number  of  unknown  nodal  displacement  components. 

Of  the  many  techniques  developed  for  solving  nonlinear  equations, 
those  based  on  the  Newton-Raphson  method  are  the  most  popular  and  widely 
known.  The  Newton-Raphson  method  was  chosen  for  use  in  this  work  because 
of  its  applicability  to  problems  involving  load  path  dependent  material 
behavior,  finite  elastic  deformations,  and  to  problems  with  both  geometric 
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and  material  nonlinearities.  The  method  is  by  no  means  perfect;  cases  exist 
where  it  diverges.  Two  examples  are  "locking"  structures  and  structures 
that  exhibit  "snap- through"  buckling,  i.e. , the  solution  process  cannot 
trace  the  unstable  portion  of  the  load-deflection  curve.  However,  it 
readily  solves  equations  arising  in  the  majority  of  practical  nonlinear 
analyses.  A brief  description  is  included  herein;  additional  details 
are  given  in  Ref.  [15]. 

> 
I 

3.4.1  The  Newton-Raphson  Method 

The  Newton-Raphson  process  forms  the  basis  of  several  proposed 
solution  strategies.  The  technique  adopted  in  this  work,  termed  incre- 
mental-iterative, is  the  most  general  and  contains  other  forms  as  special 
cases.  The  total  load,  {P^},  applied  to  the  structure  is  divided  into  a 
number  of  increments  or  load  steps  (aP^}.  The  number  and  size  of  load 
steps  chosen  for  analysis  is  usually  based  on  the  severity  of  the  nonlinear 
response  and  the  load  levels  at  which  results  (displacements,  stresses, 
and  strains)  are  desired.  The  solution  proceeds  in  a stepwise  fashion 
for  successive  load  steps  beginning  with  step  one.  Changes  of  nodal 
displacements  within  a load  step  are  computed  by  a series  of  linear 
approximations  (iterations).  The  correct  displacements  are  assumed  to 
have  been  reached  when,  for  example,  terms  of  the  residual  load  vector, 

Eq.  3.3.2,  are  sufficiently  small.  Solution  for  the  next  step  is  then 
begun. 

A major  component  of  the  Newton-Raphson  solution  process  involves 
correcting  the  nodal  displacements  to  eliminate  the  residual  loads.  For 
a given  set  of  nodal  displacements,  (Ul,  that  do  not  satisfy  the  require- 
ment of  {R}  = {0},  a correction  {dR}  is  sought  such  that 
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{R}  + {dR}  = {0} 


(3.4.1) 


A differential  change  in  the  residual  nodal  loads  is  given  by 


{dR}  = {dPj}  - {dIFj}  = [J3]  {dU} 


(3.4.2) 


where  the  Jacobian  matrix,  [J^].  contains  partial  derivatives  of  the 
function  {R}  with  respect  to  nodal  displacements  evaluated  at  the 


current  displacements.  The  Jacobian  is  conveniently  split  into  the  sum 


of  two  matrices  such  that 


{dR}  = [K]{dU}  = [ [K^]  - [K^]  ] {dU} 


(3.4.3) 


[Kj^]  is  the  "initial  load  stiffness"  and  accounts  for  changes  in 
displacement  dependent  loadings  (nonconservative  loads).  The  matrix 
[ICj.],  termed  the  "tangent  stiffness  matrix",  determines  changes  in 
internal  nodal  forces  corresponding  to  changes  in  nodal  displacements. 


Eq.  3.4.3  represents  a set  of  linear  simultaneous  equations  linking 


differential  changes  in  the  residual  load  vector  and  nodal  displace- 


ments. Substitution  of  Eq.  3.4.3  into  3.4.1  and  changing  from  differen- 


tial to  finite  increments  yields  an  expression  for  the  corrective  nodal 


displacements  that  el ini nates  the  residual  load 


{AU}  = - [K]"'{R} 


(3.4.4) 


Solution  of  the  linear  equations  is  performed  by  Gaussian  elimination 


or  some  other  triangulation  scheme  rather  than  formal  inversion  as 


indicated  above. 


A graphical  summary  of  the  various  algorithms  based  on  the  NeMton- 
Raphson  procedure  for  a single  nonlinear  equation  Is  given  In  Fig.  3.4.1. 

In  the  classical  formulation.  Illustrated  In  Fig.  3.4.1b,  [K]  Is 
regenerated  and  triangulated  before  each  Iteration  to  assure  the  best 
possible  convergence  rate.  However,  this  procedure  Is  prohibitively 
expensive  for  structures  with  many  nodes.  In  a nx)d1f1ed  version  of 
Newton- Raphson,  [K]  Is  reformed  and  triangulated  before  each  step  but 
held  constant  for  all  Iterations  within  the  step  as  shown  In  Fig.  3.4.1c. 
More  iterations  are  required  but  the  cost  of  each  iteration  is  dras- 
tically reduced.  In  the  limiting  case,  the  conventional  linear  stiff- 
ness matrix  is  used  for  all  load  steps  (Fig.  3. 4. Id).  This  modification, 
termed  the  "constant  stiffness  method",  converges  only  for  problems 
with  slightly  nonlinear  behavior,  and  even  these  may  require  an  excessive 
number  of  Iterations.  Still  other  forms  of  Newton-Raphson  update  the 
stiffness  matrix  before  specific  iterations  of  a load  step  In  an  attempt 
to  optimize  the  convergence  rate.  But  as  Is  evident,  all  the  various  forms 
are  special  cases  of  the  original  Newton-Raphson  method  and  are  easily 
Incorporated  Into  a general  solution  algorithm. 

3.4.2  Outline  of  Solution  Procedure 

The  following  is  a summary  of  the  major  computational  steps 
required  to  analyze  the  structure  for  each  load  step: 

1)  Compute  the  equivalent  nodal  loads,  corresponding 

to  the  Increment  of  applied  loads  defined  for  the  step. 

Set  {R}  = {aPj}. 

2)  If  required,  update  the  [K]  matrix  and  triangulate. 
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3)  Update  the  total  nodal  loads  currently  applied  to  the 
structure. 

4)  Solve  for  an  increment  of  the  nodal  displacements  using  the 
current  triangulated  stiffness  {aU}  = [K]"^  {R}. 

5)  Compute  increments  of  Green's  strain  and  the  new  total  2nd 
Piola-Kirchoff  stresses  for  each  nonlinear  element.  For  large 
displacement  analysis,  [B^]  is  evaluated  at  {U}  + j {aU}  to 
compute  the  correct  increment  of  strain. 

6)  Update  the  total  nodal  displacements,  element  strains,  and 
stresses  to  reflect  the  previously  computed  displacement 
increment. 

7)  Evaluate  the  internal  nodal  forces  for  all  elements  using  current 
total  displacements  and  stresses.  For  linear  elements  and 
substructures  {IF}  = [K]{U}  where  [K]  is  the  conventional  linear 
stiffness  matrix.  For  nonlinear  elements,  internal  forces  are 
computed  by  Eq.  3.2.16. 

8)  Compute  the  structure  residual  nodal  load  vector  as  {R}  = {P^}- 
{IFs). 

9)  Apply  convergence  tests  on  (R),  {U},  or  {AUl.  If  the  convergence 
test(s}  are  satisfied,  go  to  step  (1)  and  begin  processing  the 
next  step;  otherwise  go  to  step  (10). 

10)  Update  and  triangulate  the  stiffness  matrix  if  required.  Then 
return  to  step  (4)  to  begin  the  next  iteration. 

3.4.3  Convergence  Tests 

An  important  aspect  of  the  nonlinear  solution  process  involves 
selection  of  suitable  criterion  to  indicate  the  iteration  process  has 


converged  to  a sufficiently  accurate  solution.  The  Euclidean  norm 
(square  root  of  the  summed  squares)  of  the  nodal  displacement  and  residual 
load  vectors  comprise  the  most  popular  convergence  test  quantities.  Four 
common  tests  are  listed  below: 

1)  II  {R}  II  < II  {aP}||  * tolerance 

2)  Imax  {R}^  I < ||  {aP}||  * tolerance 

3)  II  {AlD^II  < II  {AU}^||  * tolerance 

4)  |max  {aU}^|  < ||  {AU}^||  * tolerance 

where  the  subscript  refers  to  the  iteration  number. 

The  first  test  simply  compares  lengths  of  the  'n*  dimensional  applied 
load  and  residual  load  vectors.  The  second  test  compliments  the  first 
by  detecting  any  highly  localized  residual  loads  that  would  be  smoothed 
over  by  the  vector  norm  computation.  These  two  tests  are  most  applicable 
to  solid  mechanics  problems  in  which  all  generalized  forces  are  of  the 
same  class.  In  frame,  plate,  and  shell  structures,  generalized  forces 
consist  of  moments,  shears,  and  membrane  forces  which  often  vary  by  several 
orders  of  magnitude  in  size.  The  first  two  convergence  tests  are  less 
appropriate  for  these  types  of  structures. 

The  third  and  fourth  tests  are  based  on  the  reduction  of  corrective 
displacement  magnitudes  as  iterations  proceed  for  a load  step.  Both  assume 
that  large  changes  in  displacement  occur  during  the  first  few  iterations 
and  that  changes  become  successively  smaller  in  subsequent  iterations. 

The  selection  of  tolerance  values  continues  to  be  a matter  of 
engineering  judgment.  The  convergence  test  and  tolerance  value  chosen 
tend  to  be  highly  problem  dependent.  Often,  a nodal  displacement  and 
residual  force  criterion  are  specified,  both  of  which  must  be  satisfied 


for  convergence.  Hand  [28]  used  a test  similar  to  the  fourth  one  above 
with  a tolerance  of  5%.  Nayak  and  Zienkiewicz  [50]  relied  on  tolerances 
ranging  from  1%  to  0.01%  for  the  first  test  with  the  total  load  norm 
substituted  for  the  incremental  load  norm.  This  latter  test  based  on  the 
total  load  vector  norm  can  be  misleading.  As  the  applied  load  on  the 
structure  increases,  the  total  load  no'*m  increases  quadratically  and  for 
a constant  tolerance  the  overall  convergence  criterion  becomes  less  severe 
especially  during  the  final  stages  of  analysis  when  nonlinear  behavior  is 
most  pronounced. 

During  this  study  numerous  structures  were  analyzed  using  the  conver- 
gence criterion  outlined  above.  The  residual  load  criterion  proved  useful 
as  a guide  for  selecting  subsequent  load  step  sizes.  It  was  also  found 
that  unduly  small  tolerances  produced  only  minor  improvements  in  the 
solution.  For  example,  differences  in  solutions  were  insignificant  for 
1%  and  0.1%  tolerances  with  the  first  convergence  test. 

3.5  Element  and  Structure  Stiffnesses 

Expressions  for  the  element  tangent  stiffness  matrix  [ICj.]  are  derived 
from  the  internal  forces 


{IF}  = 


.V 


[B]^{o}  dV 


(3.5.1) 


where  both  [B]  and  {a}  are  indirectly  functions  of  the  element  nodal 
displacements.  Differentials  of  internal  forces  are  then 


{dIF}  = 


d[B]'’^{a}  dV  + [B]^{d0}  dV 
V 


(3.5.2) 


Expansion  of  the  second  integral  using  the  material  tangent  modulus 
matrix,  Eq.  2.6.2  yields 


[ef  {do}  dV  = [ f [B^CDtICB]  dV]  {dU}  (3.5.3) 

Jv  Jv  ^ 

For  small  displacement,  linear  and  material  nonlinear  analysis,  the  [B] 
matrix  is  independent  of  nodal  displacements  (contains  only  element  nodal 
coordinates);  thus,  the  first  integral  of  Eq.  3.5.2  vanishes  leaving  the 
familiar  second  integral  as  the  element  tangent  stiffness  matrix.  For 
geometric  nonlinear  analysis  the  first  integral  must  be  included.  Un- 
fortunately, the  expansion  of  d[B]^{o}  requires  considerable  manipulation 
to  obtain  a useful  form.  Nayak  [49]  was  the  first  to  present  a matrix 
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formulation  for  a general  three  dimensional  element  and  his  notation  is 
followed  here.  A different  and  more  detailed  derivation  of  the  matrix 
formulation  is  given  in  Appendix  A.  As  shown  in  the  Appendix,  the  first 
integral  of  Eq.  3.5.2  can  always  be  written  in  the  form 

d[B]^{o}  dV  = [ f [G]^[M][6]  dV]  {dU}  (3.5.4) 

Jv  Jv 

where  [G]  contains  only  element  shape  function  derivatives  and  [M]  is  a 
symmetric  matrix  of  stresses  within  the  element.  The  above  integral 
generates  what  is  commonly  termed  the  "initial  stress"  stiffness  matrix. 
The  element  tangent  stiffness  matrix  containing  all  geometric  and  material 
nonlinear  terms  can  be  written  as 


[K^]  = [ ^[Gf[M][G]  + [Bf[D^][B]  dV]  (3.5.5) 

The  tangent  stiffness  matrix  above  is  symmetric  for  a symmetric  [Dy] 
matrix. 

The  initial  load  stiffness  matrix,  [K|^],  for  an  element  reflects 
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loadings.  Such  loads  include  pressures  on  highly  deformable  membranes 
in  which  significant  changes  of  area  magnitude  and  direction  occur.  This 
relationship  is  expressed  by 


{dP}  = [K^]  {dU} 


(3.5.6) 


Nayak  has  shown  that  [Kj^]  is  generally  a nonsymmetric  matrix  with  terms 
quite  small  compared  to  [Ky].  Details  of  the  derivation  follow  the  same 
procedure  as  used  to  derive  [Kj].  The  nonsynwetry  is  anticipated  by 
noting  that  area  transformations  between  initial  and  deformed  configura- 
tions derived  in  Chapter  2 are  nonsymmetric.  The  nonsymmetry  of  [Kj^] 
considerably  complicates  determination  of  corrective  displacements  in  the 
solution  process  by  making  the  structure  stiffness  nonsymmetric.  Non- 
symmetry  of  the  equations  essentially  doubles  the  solution  effort  required. 
The  smallness  of  terms  in  [Kj^]  combined  with  its  nonsymmetry  and  the  fact 
that  in  most  structures  loads  are  independent  of  displacement  suggest 
neglecting  [K^^]  in  computing  the  element  stiffness.  The  effect  of  any 
nonconservative  loads  is  then  accounted  for  during  computation  of  equivalent 
nodal  loads  at  the  beginning  of  each  load  step  using  the  deformed  geometry. 

In  the  Lagrangian  formulation,  once  element  stiffness  matrices  are 
computed,  assembly  of  the  structure  stiffness  matrix  follows  the  direct 
stiffness  method  as  for  linear  analysis.  This  process  is  described  in  all 
standard  texts. 


3.6  Substructuring  and  Static  Condensation 


Substructuring  and  static  condensation  are  natural  extensions  of 
the  basic  finite  element  analysis  method.  The  analyst  divides  the  structure 
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into  substructures  along  convenient  internal  boundaries  separating  logical 
components  of  the  structural  model.  Substructures  may  be  further  subdivided 
into  simple  finite  elements  or  into  other  substructures  continuing  for  as 
many  levels  as  practical.  Typically,  each  substructure  models  some  portion 
of  the  actual  structure  that  appears  repeatedly  in  the  next  higher  level 
substructure. 

A hierarchial  tree  as  shown  in  Fig.  3.6.1  provides  a graphical 
representation  for  a system  of  substructures.  The  final  structure  is 
the  root  "node"  of  the  tree;  substructures  appear  as  intermediate  nodes. 
Simple  finite  elements  must  compose  all  terminating  nodes  of  the  tree, 
i.e.,  nodes  with  no  lower  "branches".  Assembly  of  element  and  substructure 
stiffness  matrices  and  load  vectors  begins  with  the  terminating  nodes  and 
proceeds  upward  through  the  hierarchy.  Solution  for  the  nodal  displacements 
is  performed  only  for  the  highest  level  structure.  Nodal  displacements 
within  substructures  are  determined  by  a mapping  process  beginning  at  the 
highest  level  structure  and  proceeding  down  the  tree. 

In  linear  analysis,  substructuring  eliminates  the  wasteful  computation 
of  stiffness  matrices  and  equivalent  nodal  loads  for  identical  elements 
and  substructures  (those  appearing  at  more  than  one  location  in  the  tree). 
However,  all  nodes  in  lower  level  substructures  also  appear  in  higher 
level  structures.  The  final  structure  has  the  same  number  of  nodes  as  a 

t 

model  without  substructuring. 

The  number  of  nodes  in  the  highest  level  structure  can  often  be 
drastically  reduced  through  static  condensation  of  lower  level  substruc- 
tures. Consider,  for  example,  the  arbitrary  structure  shown  in  Fig.  3.6.2. 
Nodes  of  each  substructure  are  classified  as  either  interior  or  exterior. 


Interior  nodes  are  defined  as  nodes  which  do  not  directly  connect  to  nodes 
of  other  substructures.  Static  condensation  of  a substructure  produces 
a mathematically  equivalent  stiffness  matrix  and  load  vector  for  the 
exterior  nodes  .hat  contain  all  effects  of  interior  nodes.  The  actual 
substructure  is  thus  replaced  in  the  hierarchy  by  an  equivalent  substruc- 
ture with  only  exterior  nodes.  Przemieniecki  [57]  describes  the  necessary 
computations  in  detail.  The  primary  difference  between  the  classical 
formulation  given  in  the  reference  and  the  one  used  in  practice  is  that 
no  matrix  inversions  are  performed.  Instead  special  triangulation  schemes 
based  on  one  of  the  standard  decomposition  algorithms  are  employed. 

Condensed  substructures  are  introduced  into  the  hierarchy  between 
the  substructure  to  be  condensed  and  the  substructure  in  which  it  appears 
as  an  element.  Thus,  a node  in  the  tree  representing  a condensed  sub- 
structure has  only  one  lower  level  connecting  branch  as  shown  in  Fig.  3,6.1. 
When  a condensed  substructure  is  encountered  during  stiffness  or  load 
vector  assembly,  the  process  is  temporarily  suspended  while  the  condensation 
is  performed.  The  normal  assembly  process  then  resumes.  Viewed  in  this 
manner,  condensation  at  any  number  of  levels  presents  no  difficulties  in 
the  formulation.  Displacements  inside  condensed  substructures  are  obtained 
by  the  process  known  as  back-condensation  during  traversal  down  the  tree 
after  solution  of  the  highest  level  structure. 

A flexible  substructuring  and  condensation  capability  enables  the 
analyst  to  greatly  reduce  the  computer  cost  incurred  during  a nonlinear 
analysis.  Linear  regions  of  the  structure,  often  determined  by  inspection, 
are  substructured  and  condensed  leaving  only  nodes  that  interact  with 
nonlinear  substructures.  Reductions  in  computer  time  result  if  a signifi- 
cant number  of  nodes  can  be  eliminated  from  the  incremental -iterative 
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nonlinear  solution  process.  Savings  of  computer  time  occur  primarily  < 

i 

during  the  following  solution  phases:  ] 


elements  and  nonlinear  substructures  must  be  recomputed  and 


assembled.  Stiffness  matrices  of  linear  elements,  substructures. 


and  condensed  linear  substructures  computed  during  the  first 
load  step,  are  retained  for  use  throughout  the  analysis. 
Trianquiatlon  of  the  tangent  stiffness.  Most  nonlinear  analyses 
require  numerous  updates  and  tri angulations  of  the  tangent  stiffness 
of  the  highest  level  structure.  The  cost  of  analysis  Is  generally 
dominated  by  the  time  required  to  triangulate  the  equations. 
Condensation  of  linear  substructures  eliminates  many  nodes  from 
the  highest  level  structure  thereby  reducing  the  effort  required 
to  repeatedly  triangulate  the  equations.  The  time  required  to 
perform  the  initial  condensation  of  linear  substructures  is 
generally  trivial  compared  to  the  savings  achieved  In  triangulations 
of  the  nonlinear  structure  stiffness.  As  the  frequency  of 
stiffness  updates  Increases,  for  example,  near  the  ultimate  load, 
the  efficiency  of  the  solution  with  condensation  greatly  increases. 
Strain  and  stress  computations.  Since  the  behavior  of  a linear 
element  is  not  a function  of  the  strains  or  stresses,  it  is  not 
necessary  to  compute  this  data  for  linear  elements,  substructures, 
and  condensed  substructures  during  the  nonlinear  Incremental 
solution.  Displacements,  strains,  and  stresses  Inside  linear 
substructures  are  computed  only  at  the  user's  request  after  the 
nonlinear  solution  has  converged. 
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4)  Residual  Toads  calculation.  The  internal  force  vector  for 
linear  elements  and  substructures  is  computed  from  the  product 
of  the  stiffness  matrix  and  total  displacements.  This  is 
particularly  effective  for  condensed  substructures  as  only  the 
internal  force  vector  for  the  boundary  nodes  is  required. 
Residual  loads  are  zero  inside  linear  substructures  by  defini- 
tion. 
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CHAPTER  4 


REQUIREMENTS  OF  NONLINEAR  FINITE  ELEMENT  SOFTWARE 

4.1  General 

The  matrix  formulation  presented  in  the  preceeding  two  chapters 
provides  a description  of  the  numerical  procedures  for  solving  nonlinear 
finite  element  problems.  Implementation  of  this  formulation  into  a 
software  system  for  general  purpose  analysis  is  a formidable  task.  For 
example,  man-machine  communication,  structural  modeling,  analysis  restart, 
error  recovery,  computer  resource  utilization,  flexibility,  maintain- 
ability, and  portability  are  a few  of  the  difficulties  faced  by  software 
engineers.  Although  these  factors  are  not  directly  concerned  with 
formulating  or  solving  the  governing  equations,  they  generally  determine 
the  success  or  failure  of  the  software  system.  These  topics  are  a sharp 
contrast  to  the  software  problems  of  previous  generations,  when  the 
primary  concern  was  developing  new  algorithms  suited  for  digital  computers. 

In  the  present  chapter,  only  those  aspects  of  nonlinear  finite 
element  software  directly  affecting  the  structural  analyst  are  examined. 
Particular  attention  is  given  to  the  features  and  capabilities  needed 
for  successful  large  scale,  general  purpose  nonlinear  finite  element  soft- 
ware. Specific  techniques  employed  to  implement  the  type  of  software 
discussed  here  are  presented  in  Chapter  5. 

4.2  User-Program  Interface 

The  most  commonly  overlooked  aspect  of  engineering  software  involves 
communication  between  the  user  and  the  computer  program.  Input  and  output 
are  usually  viewed  by  program  developers  as  nuisances  to  be  hastily  dealt  with 
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They  fail  to  recognize  that  a user  experiences  only  three  aspects  of  a 
computer  system.  These  are  a)  data  entry,  b)  output  capability,  and 
c)  costs.  Most  users  are  willing  to  pay  extra  for  a program  that  has 
flexible  input  and  output  features. 

Traditionally,  finite  element  programs  have  been  executed  in  batch 
mode;  data  entry  has  generally  been  fixed  format.  Fixed  format  input 
is  tedious  and  prone  to  errors.  Checking  is  difficult  because  no 
descriptive  labels  appear  within  the  data.  Even  experienced  users  must 
often  refer  to  program  manuals  for  data  ordering  and  format  fields. 

Widespread  use  of  timesharing  computers  for  interactive  processing 
prompted  development  of  the  so  called  quest  answer  data  entry  mode, 
in  which  the  user  is  asked  to  respond  to  a series  of  questions  via  a 
keyboard.  The  primary  objection  to  this  type  of  input  is  that  the  user 
can  only  respond  to  the  data  requests  presented  by  the  program.  Diffi- 
culties arise  when  mistakes  are  made.  No  convenient  techniques  exist 
with  question-answer  input  to  request  that  a particular  question  be 
re-asked.  The  only  advantage  of  this  mode  of  data  entry  is  the  ease  with 
which  it  is  programmed. 

Interactive  graphics  is  becoming  a popular  means  of  data  entry  for 
describing  finite  element  models.  Generally,  a combination  of  two 
techniques  is  used.  In  the  first,  a menu  is  placed  on  the  screen  from 
which  the  analyst  may  select  one  or  more  items  via  a light  pen  or  through 
cross  hair  alignment  and  function  buttons.  Based  on  the  item  selected, 
the  user  is  shown  another  more  detailed  menu  or  he  is  asked  to  answer  a 
series  of  questions  through  a standard  keyDoard.  Although  requiring 
considerably  less  effort  on  the  analyst's  part,  these  techniques  are 
essentially  a graphical  form  of  the  typical  question-answer  input  mode. 
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The  ideal  mechanism  for  user-program  interface  in  batch,  interactive, 
and  graphic  environments  would  be  voice  connunication  via  natural  languages. 
Unfortunately,  natural  language  translators  are  in  the  early  research 
stage  and  are  many  years  away  from  practical  application.  In  the  interim, 
artificial  languages,  known  as  Problem  Oriented  Languages  (POLs),  are 
being  developed  by  software  engineers.  POLs  are  simple  languages  designed 
for  specific  engineering  disciplines.  They  consist  of  easily  remembered, 
English- like  phrases.  POLs  permit  users  to  "tell"  the  program  what  to 
do,  placing  the  user  on  the  offensive,  as  opposed  to  his  defensive 
position  in  the  question-answer  mode  of  data  entry. 

POLs  differ  from  fixed  format  and  question-answer  input  in  a number 
of  desirable  ways.  First,  commands  are  always  free  form;  words  and  data 
items  may  be  placed  anywhere  on  an  input  line.  Secondly,  most  well 
designed  POL  systems  permit  any  logical  ordering  of  commands  and  thoroughly 
check  data  for  consistency  during  the  input  phase  of  solution.  Most 
often,  corrections  can  be  made  immediately  in  an  interactive  environment 
by  simply  re-entering  the  command. 

Some  engineers  have  expressed  concern  that  translation  of  POLs  is 
inefficient  and  increases  computer  costs.  This  is  true  only  for  analysis 
of  very  small  linear  structures.  Experience  has  shown  that  for  practical 
linear  analyses  and  for  all  nonlinear  analyses,  numerical  computations 
are  the  dominating  factor  in  computer  charges  — the  cost  incurred 
translating  the  data  is  trivial. 


u 

u 


In  large  linear  analysis  systems,  batch  operation  with  fixed  format 
data  has  always  been  standard.  However,  considerable  interaction  between 
the  analyst  and  the  software  is  mandatory  for  efficient  nonlinear  analysis. 
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Engineers  must  be  more  intimately  aware  of  the  structural  response  and 
reflect  this  knowledge  in  parameters  controlling  the  solution  process. 
Blind  "number  crunching"  as  often  occurs  in  linear  analysis  is  not 
appropriate.  To  support  this  type  of  interaction,  the  software  must 
permit  a variety  of  input  modes  --  batch,  graphical,  and  interactive. 

Batch  and  graphical  input  through  a POL  is  suitable  for  specifying  the 
structural  topology  and  loads.  Multiple  analysis  restarts  via  interactive 
sessions  supported  by  a POL  are  appropriate  for  requesting  nonlinear 
analysis,  output,  and  modifying  solution  parameters. 

4.3  Structural  Modeling 

Developing  an  appropriate  finite  element  model  for  a structure  and 
specifying  that  model  to  the  computer  program  is  a time  consuming  task 
for  structural  analysts.  Finite  element  software  should  provide  features 
to  permit  adequate  modeling  of  the  nonlinear  structure  with  a minimum 
of  effort.  Several  aspects  of  this  process  are  discussed  in  the  following 
sections. 

4.3.1  Substructuring  and  Static  Condensation 

Substructuring  and  static  condensation  have  proven  effective  in 
reducing  input  data  and  computer  costs  associated  with  linear  analysis 
[4,17].  Even  greater  savings  are  possible  in  nonlinear  analysis  when 
geometric  and  material  nonlinearities  can  be  isolated  in  specific  regions 
of  a structure.  Linear  regions,  usually  determined  by  inspection,  are 
substructured  and  condensed  leaving  only  boundary  nodes  that  interact 
with  nonlinear  substructures.  The  result  is  a drastically  reduced  number 
of  nodes  present  in  the  highest  level  structure  on  which  the  incremental- 
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iterative  solution  Is  performed.  After  the  nonlinear  solution  converges 
for  the  highest  level  structure,  displacements,  strains,  and  stresses 
Inside  condensed  lower  level  linear  substructures  need  be  computed  only 
If  requested  by  the  analyst.  Substructuring  and  condensation  as  just 
described  are  demonstrated  In  two  example  problems  In  Chapter  6. 

The  bulk  of  Input  data  describing  a linear  or  nonlinear  structure 
consists  of  nodal  coordinates  and  element  Incidences  that  define  the  size, 
orientation,  and  connectivity  of  elements.  Substructuring,  combined  with 
sophisticated  mesh  generators,  reduces  the  amount  of  Input  data  for 
structures  that  contain  repeated  Identical  components.  A substructure 
is  defined  to  encompass  all  elements  and  nodes  of  a common  structural 
component.  Nodal  coordinates  In  other  occurrences  of  the  same  substruc- 
ture are  not  required  Input.  Only  the  incidences  and  orientation  of  the 
substructure  as  It  appears  in  the  higher  level  structure  are  necessary 
[34,  35].  Rotation  angles  are  a convenient  means  of  describing  a 
substructure's  orientation. 

Similarly,  various  patterns  of  loads  applied  to  nodes  and  elements 
of  a substructure  are  defined  only  once.  Sets  of  effective  nodal  loads 
may  be  selectively  applied  to  each  occurrence  of  the  substructure  In  a 
higher  level  structure  to  model  the  real  applied  loads.  The  orientation 
of  loads  on  each  occurrence  normally  follows  that  of  the  substructure; 
however,  an  engineer  can  potentially  specify  load  orientation  Independently 
of  substructure  orientation.  A simple  example  Involves  gravity  type 
loads  In  which  the  direction  remains  constant  even  though  the  sub- 
structure orientation  is  altered. 
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For  maximum  effectiveness,  substructuring  and  condensation  must  be 
an  integral  feature  of  the  structural  definition  process  from  the  user's 
viewpoint.  Once  defined,  a substructure  used  in  another  structure  should 
appear  as  any  other  finite  element  to  the  analyst  during  input.  Computa- 
tional details,  such  as  mapping  nodal  degrees  of  freedom  between  sub- 
structures, should  be  handled  completely  by  the  software.  The  analyst 
should  not  be  required  to  specify  the  order  and  kind  of  operations  necessary 
to  effect  multiple  level  condensations  or  to  specify  structural  data 
that  must  be  recomputed  during  a nonlinear  analysis.  No  currently  available 
nonlinear  analysis  system  provides  a completely  automated  substructuring 
and  condensation  capability.  Only  a few  linear  systems  offer  any  type  of 
substructuring.  Implementation  of  these  features  greatly  increases 
program  complexity  and  development  time  if  multi-level  condensations  with 
loads  are  permitted.  In  spite  of  the  obvious  complexities  introduced 
into  the  software  to  support  substructuring  and  condensation,  these 
features  appear  essential  for  efficient  solution  of  large  nonlinear 
structures. 

4.3.2  Loads 

Loading  conditions  defined  for  linear  analysis  generally  consist  of 
nodal  forces,  element  pressures,  temperature  gradients,  etc.  that  corres- 
pond to  some  real  loading  situation.  In  nonlinear  analysis,  the  concept 
of  a loading  condition  is  extended  to  imply  a loading  "pattern"  which 
defines  a spatial  distribution  of  loads  on  a structure.  A single  nonlinear 
load  step  may  consist  of  any  number  of  loading  patterns  combined  with 
scalar  multipliers  to  form  the  real  load  increment.  To  incorporate  time 
dependent  effects,  such  as  creep  and  shrinkage,  the  analyst  must  be  able  to 
associate  a time  parameter  with  the  definition  of  a nonlinear  load  step. 
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Multiple  load  patterns  in  a load  step  are  essential.  For  instance, 
analysts  may  want  to  simulate  a construction  sequence  by  applying  dead 
loads  followed  by  live  loads  in  different  steps.  Similarly,  certain 
loadings  may  be  simpler  to  describe  when  broken  into  components  and 
applied  as  different  loading  patterns.  Systems  that  permit  only  one 
pattern  in  a load  step  limit  the  analysis  to  proportional  loading. 

Loading  definition  becomes  more  complex  with  the  introduction  of 
multi-level  substructuring  and  condensation.  The  following  technique 
developed  during  this  study  appears  to  be  quite  flexible. 

As  lower  level  substructures  in  the  hierarchy  are  defined,  the 
declaration  of  loading  patterns  is  performed  as  for  linear  loading 
conditions.  In  higher  level  substructures,  loading  patterns  are  defined 
in  terms  of  patterns  on  lower  level  substructures.  (Lower  level  sub- 
structures are  treated  as  elements  in  higher  level  structures  during 
structural  definition).  Finally,  any  number  of  loading  patterns  may 
exist  on  the  highest  level  structure.  A special  type  of  loading  condition, 
designated  as  nonlinear,  is  declared  only  for  the  highest  level  structure. 

A nonlinear  loading  condition  specifies  individual  step  load  increments 
consisting  of  multiples  of  pattern  loads  on  the  highest  level  structure 
with  any  additional  time  parameters. 

During  solution  for  a load  step,  the  equivalent  nodal  loads  for  each 
pattern  specified  in  the  step  are  computed  and  combined  with  the  specified 
multipliers.  For  loading  patterns  that  appear  in  more  than  one  step. 


the  equivalent  nodal  loads  need  to  be  computed  only  once  and  saved  between 


load  steps.  However,  equivalent  nodal  forces  for  loading  patterns  on 


nonconservative  effects.  The  software  should  automatically  determine  the 
appropriate  sequence  of  operations  (including  condensation  of  loads)  to 
effect  the  proper  solution  without  intervention  by  the  analyst. 


4.3.3  Constraints 


The  tangent  stiffness  matrix  of  a structure  contains  rigid  body 
motions  that  prevent  a unique  solution  for  displacements.  A sufficient 
number  of  constraint  equations  to  remove  all  possible  rigid  body  motions 
are  necessary  before  the  equilibrium  equations  are  rendered  non-singular 
Any  number  of  additional  consistent  and  nonredundant  constraints  are 
permitted  to  model  the  physical  boundary  conditions  of  the  structure. 

Two  classes  of  constraints  are  necessary  for  general  analysis.  The 
first  class,  termed  absolute  constraints,  are  of  the  form 


U.  = constant 


and  force  the  i degree  of  freedom  displacement  to  equal  the  constant 


The  second  class,  termed  relative  constraints,  force  a linear  relation 


ship  between  displacements  of  two  or  more  degrees  of  freedom.  Such 
constraints  are  expressed  by 


Zci^.U^.  = constant 


Relative  constraints  are  useful  for  connecting  various  types  of  finite 


elements,  for  instance,  beam  and  shell  elements 


During  a nonlinear  incremental  solution,  constraints  also  are  incre 
mental.  Thus,  absolute  and  relative  constraints  for  nonlinear  analysis 


are  expressed  by 


aU<  = constant 


Za,AU.  = constant 


0 
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A special  case  occurs  when  the  constant  is  always  zero.  Degrees  of  free- 
dom with  absolute  constraints  then  always  have  zero  displacements.  With 
relative  constraints,  the  relationship  between  displacements  remains 
unaltered  throughout  the  analysis. 

The  interpretation  of  constraints  with  non-zero  constants  must  be 
incremental.  During  the  first  Iteration  of  a load  step,  corresponding 
to  an  Increment  of  applied  loads.  Incremental  displacements  are  forced 
to  satisfy  the  constraint  equations.  If  further  Iterations  are  necessary 
for  a step  to  correct  the  linearized  displacements,  the  constraint 
equation  constants  should  be  made  temporarily  zero.  Otherwise,  at 
completion  of  the  Iterations  the  sum  of  displacement  changes  occurring 
over  the  step  would  not  satisfy  the  specified  constraint  equations  for 
the  step. 

4.4.  Solution  Algoritfmis  and  Convergence  Criterion 

The  Newton-Raphson  procedure  is  at  present  the  most  widely  applicable 
method  for  solving  the  nonlinear  finite  element  equilibrium  equations. 

The  procedure  has  a very  simple  physical  Interpretation  and  converges 
relatively  fast  for  most  problems.  Numerical  analysts  have  not  found 
the  method  particularly  appealing  because  a good  initial  estimate  of  the 
solution  Is  required.  However,  this  Is  seldom  a problem  in  structural 
applications  as  the  previously  converged  nonlinear  solutions  provide  a 
good  approximation  for  subsequent  load  steps. 

Several  forms  of  the  Newton-Raphson  method  have  been  developed;  all 
are  essentially  the  same  algorithm  as  discussed  In  Chapter  3.  A general 
system  must  provide  the  analyst  with  at  least  the  options  shown  in 
Chapter  3.  The  algorithm  Is  controlled  by  essentially  two  parameters: 
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a)  The  frequency  of  tangent  stiffness  updates; 

b)  The  number  of  iterations  per  step  to  improve  the  solution 
and  the  criterion  for  terminating  the  iterative  process. 

The  analysis  of  some  large  structures  may  necessitate  the  use  of 
condensed  nonlinear  substructures  to  overcome  software  size  limitations. 
Extension  of  the  iterative  technique  and  convergence  criterion  requires 
additional  consideration  because  both  the  nonlinear  substructures  and 
the  highest  level  structure  have  distinct  residual  loads  and  corrective 
displacements  (no  corrections  are  required  in  linear  substructures). 
Residual  loads  for  individual  nonlinear  substructures  are  also  reflected 
in  residual  loads  of  higher  level  substructures  as  a result  of  the 
condensation  process.  To  assure  convergence  of  the  nonlinear  solution 
in  all  parts  of  the  structure,  the  Newton-Raphson  convergence  tests 
should  be  performed  on  all  condensed  nonlinear  substructures  as  well  as 
the  highest  level  structure.  The  solution  is  said  to  converge  when  the 
tests  are  satisfied  in  all  substructures. 

4.5  Analysis  Restart  and  Error  Recovery 

A general  restart  capability  is  mandatory  for  all  forms  of  nonlinear 
analysis.  The  analyst  is  seldom  able  to  specify,  a priori,  all  loading 
increments,  convergence  parameters,  output  requests,  etc.  Users  typically 
will  discontinue  the  analysis  after  several  load  steps,  examine  struc- 
tural behavior  for  the  current  load  level,  modify  data  for  subsequent 
steps,  and  then  resume  analysis  in  another  computer  run.  This  process 
may  be  repeated  numerous  times  for  complicated  structures. 

Associated  with  restart  capability  is  the  ability  to  recover  from 
various  errors  during  solution.  Suppose,  for  example,  that  the  structure 
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suddenly  becomes  unstable  and  the  solution  diverges  after  several  load 
Increments.  If  large  load  Increments  are  used,  the  critical  load  of  the 
structure  Is  not  well  defined.  In  this  situation,  the  analyst  redefines 
the  last  load  step,  declares  additional  load  steps,  and  requests  a resump- 
tion of  the  analysis  from  the  last  converged  step.  To  accomplish  this,  a 
system  need  only  retain  the  appropriate  data  to  resume  at  any  step  based 
on  the  type  of  non! Inearl  ties  present  In  the  structure.  Retaining  all  ■ 

computed  data  for  each  step  wastes  space  on  secondary  storage  devices. 

Analysts  should  be  permitted  to  request  destruction  of  data  for  steps  no  ] 

longer  needed,  thereby  releasing  the  space  for  reuse.  | 

1 

Certain  parts  of  a structural  description,  for  example.  Incidences 
and  coordinates,  cannot  be  changed  during  a nonlinear  solution.  However, 
users  should  be  able  to  easily  modify  those  aspects  of  the  structural 
description  and  solution  parameters  listed  below: 

a)  Loading  definitions  — new  loading  patterns  and  load  step 
Increments: 

b)  Convergence  parameters  — this  Includes  both  the  types  of 
convergence  tests  and  the  tolerances; 

c)  Solution  algorithm  — the  solution  algorithm  Is  controlled  by 
the  frequency  of  stiffness  updating  and  the  number  of  equilibrium 
Iterations  permitted; 

d)  Constraints  — Incremental  constraints  can  be  modified  between 
load  Increments: 

e)  Solution  traces  — Simple  traces  of  solution  computations  permit 
the  analyst  to  determine  If  the  analysis  Is  proceeding  satis- 
factorily, or  If  changes  will  be  needed  before  the  next  load  step. 


4.6.  Element  and  Material  Model  Libraries 

During  the  lifetime  of  a finite  element  program,  system  functions, 
including  input  data  translation,  forming  and  solving  the  equilibrium 
equations,  etc.,  remain  relatively  unchanged.  Occasionally  new  features 
are  added  and  errors  corrected  but  the  basic  concepts  of  equation  solving, 
data  management,  etc.  do  not  appreciably  change.  In  contrast,  element 
formulations  and  models  of  nonlinear  material  behavior  are  still  primary 
areas  of  current  research  and  thus  change  quite  frequently.  The  library 
of  a general  system  will  be  constantly  modified  by  users  in  a research 
environment.  A measure  of  the  versatility  of  finite  element  software  is 
the  relative  ease  with  which  users  can  install  and  test  new  elements  and 
material  models. 

The  following  sections  briefly  summarize  the  function  of  element 
and  material  model  libraries  and  mechanisms  for  interfacing  library 
routines  with  the  system. 

4.6.1  Finite  Element  Libraries 

An  important  feature  of  the  finite  element  method  is  that  all 
elements,  regardless  of  their  complexity,  perform  essentially  the  same 
function  from  a systems  viewpoint.  During  nonlinear  analysis,  element 
dependencies  are  limited  to  the  following  specific  functions: 

a)  Generate  the  tangent  stiffness  matrix  consistent  with  the 
Lagrangian  formulation  given  the  element  nodal  displacements, 
stresses,  coordinates,  etc. 

b)  Generate  an  increment  of  Green  strain  given  the  previous  total 
displacements  and  the  increment  of  element  nodal  displacements. 
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In  finite  deformation  problems,  element  strain  routines  must 
also  generate  a matrix  of  geometric  transformation  data  that 
permits  material  models  and  output  routines  to  transform  between 
various  definitions  of  strain  and  stress. 

c)  For  linear  or  geometrically  nonlinear  elements,  the  element 
stress  routine  computes  total  stresses  given  the  previous 
stresses  and  the  Green  strain  Increment.  In  materially  nonlinear 
analysis,  element  stress  routines  typically  perform  no  computa- 
tions. Exceptions  are  some  beam,  plate,  and  shell  elements, 

for  example,  that  Integrate  stress  conq)onents  to  form  stress 
resultants. 

d)  Pressures,  temperature  gradients,  and  other  applied  element  loads 
are  converted  to  equivalent  element  nodal  loads  by  element  load 
routines.  Effects  of  nonconservative  loads.  If  present,  are 
accounted  for  In  these  routines. 

e)  During  residual  load  computation,  element  routines  generate 
the  Internal  nodal  forces  corresponding  to  the  current  element 
stresses  and  deformed  geometry. 

f)  Some  elements  may  generate  special  output  data,  such  as  principal 
stresses,  or  alter  the  number  of  points  for  output.  This  data 

Is  not  Incremental  In  nature  and  therefore  cannot  be  computed 
during  solution.  These  values  must  be  computed  by  element 
routines  just  prior  to  output. 

4.6.2  Material  Model  Libraries 

Material  models  are  necessary  to  Idealize  the  behavior  of  real 
materials  under  complex  states  of  stress  In  the  nonlinear  range.  Incremental 


• f-Tn- 


66 


deformation  theories  capable  of  modeling  a broad  range  of  behaviors  are 
appropriate  for  general  analysis.  Model  dependencies  can  be  limited  to 
the  following  operations: 

a)  Initialize  history  dependent  quantities  that  describe  subsequent 
material  behavior.  Form  the  elastic  stress-strain  relationship 
constitutive  matrix. 

b)  Compute  the  new  total  stresses  at  a point  within  an  element 
given  the  previous  stresses,  strain,  history  parameters,  etc. 
and  the  Increment  of  Green  strain.  Effects  of  creep,  shrinkage, 
temperature,  and  time  may  be  accounted  for  during  these  calcula- 
tions. 

c)  Compute  an  Instantaneous  tangent  modulus  matrix  for  re-computing 
the  element  tangent  stiffness  matrix. 


4.6.3  System-Library  Interface 

The  main  factor  determining  the  ease  with  which  users  may  modify 
the  library  of  elements  and  material  models  Is  the  Interface  between 
the  system  and  the  library  routines.  The  fundamental  requirement  of  the 
system- library  Interface  Is  standardization.  All  Interfaces  should  be 
Identical  regardless  of  the  element  or  model  complexity.  Data  made 
available  to  library  routines  should  come  directly  through  the  system 
Interface  and  not  Indirectly  through  FORTRAN  COMMON  areas  from  other 
library  routines.  The  latter  method  leads  to  confusion  for  users  not 
familiar  with  all  other  parts  of  the  system.  No  elements  or  models 
should  receive  special  consideration.  For  simple  elements  and  models 
this  leads  to  overkill  In  that  more  data  than  required  Is  made  available 
to  the  element  and  material  model  routines. 
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By  standardizing  the  Interface  and  forcing  all  data  to  pass  through 
this  Interface,  element  and  material  models  are  Isolated  completely  from 
the  system.  This  permits  the  system  to  perform  additional  functions 
coninon  to  all  elements  and  models.  Memory  allocation,  translation  of 
Input  data,  printing  of  results,  compatibility  checking,  etc.  can  all 
be  handled  by  the  system.  Moving  these  tasks  to  the  system  level  frees 
the  developer  to  pursue  his  real  Interest  — performance  testing  of  a new 
element  or  material  model. 
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CHAPTER  5 


IMPLEMENTATION  OF  FINITE  ELEMENT  SOFTWARE 


5.1  General 

Finite  element  software  is  now  recognized  by  the  engineering  pro- 
fession as  a major  link  in  the  transfer  of  technology  between  researchers 
and  practicing  engineers  [47,65,66].  Frequently,  software  represents 
the  only  practical  result  of  a research  project  since,  unlike  the 
associated  technical  publications,  it  is  already  in  a directly  usable 
form. 


Planning,  design,  and  implementation  of  finite  element  software  to 
achieve  all  the  requirements  outlined  in  the  previous  chapter  constituted 
a major  part  of  this  research.  Currently  available  nonlinear  analysis 
programs  achieve  only  a few  of  the  requirements,  primarily  because  they 
were  developed  by  finite  element  researchers  interested  in  the  formulation 
of  numerical  algorithms,  not  in  software  engineering.  As  the  size  of 
structures  increases,  inadequacies  of  these  programs  become  apparent,  as 
does  the  need  for  careful  planning  and  design  of  the  software.  The  per- 
centage of  development  effort  devoted  to  numerical  algorithms  is  relatively 
insignificant  compared  to  that  required  for  design  of  an  internal  system 
and  data  structure  capable  of  efficiently  processing  the  enormous  volumes 
of  structural  data. 

This  chapter  presents  various  design  and  implementation  aspects  of 
FINITE,  a general  system  for  linear  and  nonlinear  structural  analysis. 
FINITE  represents  a first  attempt  to  synthesize  both  the  analytical 
capability  demanded  by  sophisticated  analysts  and  modern  software  design 
technology  into  a comprehensive  nonlinear  analysis  system. 


i 

i 

L, 


Li 


;1 


f 

I 


r1 


ri 


LI 


The  present  chapter  Is  divided  Into  five  sections  for  presenta- 
tion. In  the  first,  two  major  approaches  of  software  development  are 
delineated.  The  second  section  presents  the  fundamental  concepts  embodied 
In  the  POLO  supervisor  with  examples  Illustrating  their  application  to 
nonlinear  analysis.  A brief  overview  of  FINITE,  given  In  the  third 
section,  emphasizes  the  distinguishing  characteristics  of  the  system's 
Internal  organization.  The  fourth  section  considers  nonlinear  material 
modeling  In  FINITE  from  the  user  and  new  model  developer  points  of  view. 
The  final  section  examines  two  of  the  primary  nonlinear  processing  modules 
appended  to  the  linear  version  of  FINITE  to  produce  the  nonlinear  version. 

No  attempt  has  been  made  here  to  fully  describe  details  of  the 
FINITE  system  (the  current  nonlinear  version  consists  of  more  than 
100,000  FORTRAN  statements).  Rather,  the  fundamental  concepts  of  Its 
design  are  shown  to  Illustrate  the  primary  differences  between  FINITE 
and  other  existing  nonlinear  finite  element  codes.  Chapter  6 demon- 
strates the  problem  solving  capability  of  FINITE  through  example  analyses. 
Details  of  system  usage  can  be  found  In  the  FINITE  user's  manual. 

5.2  Structure  of  Finite  Element  Software 

Interest  In  various  approaches  for  developing  finite  element 
software  arose  as  engineers  began  to  examine  the  Internal  structure  of 
computer  programs  and  associated  data.  As  a result  of  software  develop- 
ment efforts  over  the  past  decade,  two  distinctly  different  approaches 
have  emerged.  In  the  first,  only  standard  features  of  an  algorithmic 
language,  such  as  FORTRAN,  are  employed  In  the  actual  code.  Such  an 
approach  places  severe  limitations  on  the  types  of  program  and  data 
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structures  available  to  developers.  All  logical  operations  on  data  must 
be  repeatedly  programmed  by  system  developers.  The  vast  majority  of 
existing  finite  element  software  falls  within  this  category.  The 
second  approach  is  distinguished  by  the  presence  of  a comprehensive 
engineering  data  base  and  memory  management  system  that  relieves  developers 
from  programming  these  tasks.  This  approach  has  been  termed  "integrated" 
software  In  the  literature.  Characteristics  of  both  approaches  are 
discussed  In  the  following  sections. 

5.2.1  Systems  Without  a Data  Base  Manager 

Existing  nonlinear  finite  element  programs  have  a structure  similar 
to  that  shown  In  Fig.  5.2.1.  Problem  data  Is  divided  among  any  number 
of  sequentially  accessed  files  each  of  which  has  a trivial  structure.  A 
large  number  of  files  Is  necessary  (10  or  more  files  is  not  uncommon)  to 
represent  data  associated  with  a nonlinear  analysis. 

The  computer  code  consists  of  FORTRAN  subroutines  grouped  together 
into  processing  modules.  One  module  contains  the  logic  to  direct  execu- 
tion of  those  at  a lower  level.  Data  is  transmitted  between  modules  in 
two  ways.  Simple  control  parameters  are  passed  through  standard  features 
of  the  FORTRAN  language.  The  majority  of  data  Is  implicitly  passed 
through  the  sequential  files.  Each  module  reads  data  into  memory  from 
appropriate  files,  performs  the  computations,  then  re-writes  the  data  to 
the  sequential  file  prior  to  Initiating  the  subsequent  module.  Various 
schemes  that  maintain  certain  data  In  memory  at  all  times  have  also  been 
used.  The  SAP  [7,71]  family  of  codes  employs  a combination  of  sequential 
files  and  fixed  data  In  memory. 
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Inadequacies  of  this  approach  becotne  apparent  when  one  tries  to 
modify  or  expand  the  program  capabilities,  or  to  solve  large  problems. 

These  programs  become  I/O  bound  (Inordinate  amounts  of  data  transferred 
between  memory  and  files)  for  large  problems  because  each  module  performs 
Its  own  programmed  read/write  operations  Irrespective  of  the  memory 
space  available.  These  programs  cannot  utilize  additional  memory  If 
made  available,  and,  because  data  Is  tied  d1r:}ct1y  to  physical  locations 
In  a sequential  file,  changes  In  the  data  structure  usually  necessitate 
re-writing  much  of  the  program.  The  definition  of  more  files  circumvents 
this  difficulty  at  the  expense  of  further  complicating  subsequent  program 
modifications. 

The  simple  data  structures  permitted  on  sequential  files  place 
severe  limits  on  the  level  of  sophistication  achievable  In  these  programs. 
Consider,  for  example.  Implementing  a multi-level  substructuring  capability 
with  sequential  files.  The  many  structural  topologies,  geometries, 
loads,  computed  results,  etc.  would  require  an  excessive  number  of  data 
files.  With  this  approach  to  software  development.  Implementation  of 
such  features  Is  both  unrealistic  and  Impractical. 

5.2.2  Systems  With  a Data  Base  Manager 

A more  realistic  finite  element  system  Is  Illustrated  In  Fig.  5.2.2. 

A number  of  processing  modules  are  developed  that  operate  on  a logical 
data  space.  A data  management  system.  Independent  of  all  processing 
modules,  resolves  logical  data  requests  Issued  by  the  modules  Into  physical 
locations  within  a file.  Thus,  the  mapping  of  the  logical  data  space 
onto  an  actual  secondary  storage  device  Is  transparent  to  the  processing 
modules.  Each  module  may  access  any  data  within  the  data  space.  Inter- 
facing between  modules  Is  accomplished  easily  through  the  logical  data 
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•structure.  Ideally,  modifications  of  the  logical  data  structure  as  the 
system  capability  increases  should  not  require  changes  in  existing  modules. 
A memory  management  system  that  interfaces  with  the  data  manager  is 
necessary  to  map  data  from  the  physical  file  into  memory  for  use  by  the 
processing  modules. 

A much  simplified  version  of  this  ideal  scheme  has  been  implemented 
in  a number  of  large  linear  analysis  systems,  the  most  notable  of  which 
are  ASKA  and  NASTRAN.  In  these  systems,  references  to  data  management 
routines  are  imbedded  within  the  FORTRAN  code  of  the  individual  processing 
modules.  Most  references  simply  extract  a specified  matrix  from  a file. 
Hierarchial  data  structures,  required  to  implement  many  of  the  features 
described  previously,  are  difficult  to  maintain  in  this  manner.  Engineers 
developing  new  features  must  write  well  planned  sequences  of  references 
to  data  management  routines  to  interact  with  a hierarchial  data  structure. 
Additional  complications  occur  when  defining  and  maintaining  the  logical 
description  of  the  data  space  through  references  to  special  data  management 
routines. 

The  POLO-FINITE  system  more  closely  approaches  the  system  structure 
shown  in  Fig.  5.2.2.  The  POLO  supervisor  supports  the  logical  data 
space  concept  and  minimizes  programming  of  data  management  tasks  by  pro- 
viding developers  with  high  level  languages  and  compilers  to  do  the  work 
automatically.  The  FINITE  system  has  a modular  structure  centered  around 
a hierarchial,  logical  data  space.  This  philosophy  of  design  permits 
the  system  capability  to  evolve  naturally  without  major  revisions  in 
existing  modules.  A combined  interactive-batch  processing  capability  with 
POL  input  makes  the  system  flexible  and  well-suited  for  nonlinear  analysis. 
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5.3  POLO  — An  Engineering  Software  Support  System 

POLO  was  a prominant  factor  in  the  design  and  implementation  of  the 
nonlinear  capabilities  of  the  FINITE  system.  A brief  description  of 
polo's  capabilities  pertaining  to  this  study  is  presented  here.  Complete 
functional  details  of  POLO  may  be  found  elsewhere  [37]. 

POLO  provides  software  engineers  with  an  enhanced  programming 
environment  in  which  to  develop  engineering  application  software  (POLO 
subsystems).  POLO  itself  does  not  solve  engineering  problems;  rather  it 
supports  functions  common  to  all  engineering  application  areas.  These 
functions  include: 

1)  Problem  oriented  language  translation; 

2)  Data  structure  definition; 

3)  Data  and  memory  management  during  execution; 

4)  Integration  of  application  subsystems. 

POLO  supports  two  higher  level  languages,  actually  POLs,  termed 
F and  G.  Subsystem  developers  describe  the  logical  structure  of  data 
with  the  F language.  These  symbolic  data  definitions  (File  Definitions 
in  POLO  terminology)  are  converted  to  internal  form  and  stored  in  a POLO 
data  base  independently  of  application  subsystems.  After  this  process  is 
completed,  subsystem  developers  write  drivers  for  processing  modules  in 
the  POLO  host  language  G.  Statements  in  the  G language  form  "grammars" 
similar  in  function  to  a well  structured  FORTRAN  main  program.  Once 
the  grammar  is  written,  the  subsystem  is  completed  by  writing  standard 
FORTRAN  subprograms  referred  to  by  the  grammar. 

The  primary  purpose  of  a grammar  is  to  direct  execution  of  the 
supporting  subprograms  and  to  pass  them  data  from  a POLO  data  base,  i.e.. 
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the  grannar  provides  a logical  interface  between  the  data  space  and  the 
FORTRAN  subprograms.  Data  required  by  subprograms  is  referenced  in  the 
grammar  using  symbolic  names  provided  in  the  file  definition.  Instruc- 
tions for  locating  the  data  (data  management  requests)  are  generated 
automatically  by  POLO.  The  advantage  of  this  approach  is  that  subsystem 
developers  need  not  be  concerned  with  the  physical  mapping  of  data  in 
the  data  bases.  Grammars  are  also  used  to  translate  POLs  and  to  initiate 
other  POLO  subsystems.  A special  POLO  subsystem  (compiler)  converts 
the  grammars  to  an  internal  form,  termed  object  grammars,  and  saves  them 
in  a system  data  base. 

During  execution,  the  primary  POLO  interpreter  executes  the  object 
form  of  the  grammar.  When  instructed  to  invoke  a subprogram  referenced 
in  the  grammar,  POLO  routines  make  the  requested  data  physically  present 
in  memory  (if  it  is  not  already  present),  then  pass  the  data  to  applica- 
tion subprograms  as  standard  FORTRAN  arguments.  Data  is  made  present  in 
memory  by  a separate  data  management  interpreter  coupled  with  a demand 
paging  dynamic  memory  allocator.  Application  subprograms  are  not  aware 
of  data  and  memory  management  functions  as  all  operations  of  making  data 
present  are  done  prior  to  calling  them. 

In  summary,  POLO  subsystems  consist  of  three  major  components: 

1)  Symbolic  file  definitions, 

2)  Grammars, 

3)  Subsystem  FORTRAN  subprograms. 

The  following  sections  describe  the  first  two  of  these  components 
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and  illustrate  particular  applications  to  nonlinear  analysis.  Then  a 
description  of  POLO'S  run-time  configuration  is  presented. 


5.3.1  File  Definitions 

Most  data  for  finite  element  analysis  is  naturally  structured  in 
a hierarchial  manner.  Consider  a simplified  portion  of  the  FINITE  data 
structure  shown  in  Fig.  5.3.1.  The  highest  level  table,  called  STRUCTURES, 
contains  information  about  structures  and  individual  finite  elements. 

Each  row  describes  an  attribute  for  an  element  or  structure  in  a particular 
column.  Row  3 points  to  a lower  level  table,  LOADS,  that  contains  a 
description  of  the  loads,  in  addition  to  results  for  the  individual  steps 
of  a nonlinear  analysis.  Each  column  of  the  LOADS  table  contains  results 
for  a single  load  step.  Load  tables  exist  belov.'  only  these  cclunns  of 
STRUCTURES  that  contain  structures. 

One  set  of  results  examined  here  for  illustration  is  the  geometric 
transformation  matrices  between  deformed  and  undeformed  element  configura- 
tions for  large  displacement  analysis.  A matrix  of  geometric  data  is 
required  for  each  strain  point  of  every  element  that  permits  large  dis- 
placements. Row  EPSTRANSFORM  of  the  LOADS  table  points  to  a lower  level 
table  containing  pointers  to  vectors  of  transformation  data  for  many 
elements.  The  vector  number  and  position  within  the  vector  of  the  matrix 
for  the  first  strain  point  of  an  element  is  sufficient  to  locate  all 
matrices  for  the  element  (matrices  are  stored  contiguously  within  a data 
vector).  Row  EPSVECPOS  of  the  STRUCTURES  table  contains  the  vector  number 
and  position  for  the  element's  data. 

The  logical  data  structure  described  above  is  a combination  of 
« 

inverted  lists  and  hierarchial  tables.  Similar  constructs  are  used 
throughout  FINITE  to  permit  random  access  to  nonlinear  element  data 
independently  of  the  structure  in  which  the  elements  are  used.  Somewhat 
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different  data  structures  in  FINITE  permit  rapid  access  to  data  for  all 
elements  of  a specific  structure. 

The  data  structure  described  in  this  example  demonstrates  the  funda- 
mental kinds  of  data  tables  supported  by  POLO.  STRUCTURES,  LOADS,  RESULTS, 
and  EPSTRANSFORM  are  called  labelled  tables.  Each  row  may  contain 
different  types  of  data  (alpha,  real,  pointers,  etc.).  Such  tables  may 
consist  of  a single  block  of  columns  or  multiple  blocks  of  columns.  New 
blocks  of  columns  are  made  available  on  demand;  however,  grammar  references 
are  unaware  of  the  internal  blocking.  The  ETRNVEC  table  is  the  least 
flexible  table  type  available.  All  data  is  of  a single  mode  and  the  table 
is  always  contiguous  in  both  memory  and  the  data  base. 

Any  number  of  hierarchies  similar  to  the  one  shown  in  Fig.  5.3.1 
can  be  defined  within  a POLO  data  base.  The  data  structure  shown  in  the 
figure  is  defined  to  POLO  through  the  language  F with  the  following 
commands. 


FILE  DEFINITION  FINITE 
TABLE  STRUCTURES  LABELLED  GROUPING  25 
SNAME  ALPHA  2 
EPSVECPOS  INTEGER 
LOADS  LABELLED  GROUPING  5 
NAME  ALPHA  2 
RESULTS  LABELLED  1 
DISPLACEMENTS  LABELLED  GROUPING  11 

EPSTRANSFORM  LABELLED  GROUPING  11 
N INTEGER 

ETRNVEC  VECTOR  INTEGER  N 
END  OF  TABLE 
END  OF  TABLE 
END  OF  TABLE 


end'of  file  definition 
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5.3.2  Gratranars 

The  primary  function  of  a grammar  is  to  initiate  application  subprograms 
passing  them  data  from  the  data  base.  A single  reference,  no  matter  how 
deep  in  a hierarchy,  is  passed  as  one  argument  to  the  application  subprogram. 
For  example,  in  Fig.  5.3.1,  assume  the  structure  column  number  SCOL,  loading 
column  number  LCOL,  the  vector  number  for  element  transformations  VECNO, 
and  position  POS  are  all  known.  Data  for  the  element  can  be  passed  to 
the  application  subprogram  (with  symbolic  name  TRANSFORM)  via  the  following 
statement  in  the  gratmar 
USE  FINITE 
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♦EXECUTE  TRANSFORM(  STRUCTURES(LOADS,SCOL, RESULTS, LCOL, 

EPSTRANSF0RM,1,ETRNVEC,VECN0,P0S)  ) 

The  USE  command  designates  the  file  definition  and  appears  only  once  in 

the  grammar.  The  data  reference  above  causes  the  POS  entry  of  the  data 

vector  to  be  passed  to  the  subprogram  as  an  argument.  The  first  few 

lines  of  the  subprogram  might  be 

SUBROUTINE  TRNSFM(  TVECTR  ) 

DIMENSION  TVECTR(l) 

• 

RETURN 

END 

Traversals  through  the  data  hierarchy  to  resolve  a logical  data 
request  are  generally  slow.  Repeated  references  like  the  one  above  must 
be  avoided  to  attain  efficiency.  POLO  provides  additional  access  mechan- 
isms that  considerably  speed  up  execution.  Artificial  "base"  tables  can 
be  defined  at  intermediate  levels;  subsequent  lower  references  may  then 


r ' 


r 

i 


I 


I 


78 

begin  at  the  intermediate  level.  Subsystem  developers  can  manipulate 
the  actual  "pointers"  and  request  POLO  to  perform  direct  accesses  to  data 
without  the  overhead  associated  with  a traversal.  The  payoff  is  increased 
execution  speed  but  with  some  increased  programming  effort. 

5.3.3  Run-Time  Configuration 

During  execution,  POLO  and  the  various  FINITE  subsystems  appear  as 
a single  executable  FORTRAN  program  as  shown  in  Fig.  5.3.2.  The  POLO 
executive  interpreter  is  the  highest  level  driver  in  the  system.  It 
requests  the  token  scanner  to  read  POL  commands  from  various  input  devices, 
convert  free-form  data  to  fixed-format  data,  and  place  the  results  in  a 
FORTRAN  COMMON  area  accessible  by  both  POLO  system  routines  and  FINITE 
subsystem  routines.  The  top  portion  of  this  vector  is  reserved  for 
scanner  output  and  a static  area  in  which  subsystem  routines  may  place 
variables.  The  POLO  dynamic  memory  allocator  uses  the  remaining  portion 
of  the  COMMON  area  to  maintain  requested  data  in  memory. 

As  shown  in  the  figure,  application  subprograms  are  totally  isolated 
from  the  system  functions  necessary  to  locate  data  and  make  it  present 
in  memory.  Data  modified  by  the  subprograms  will  eventually  be  forced 
into  the  data  bases  either  during  execution  as  the  dynamic  allocator 
seeks  to  make  room  for  other  data  or  at  the  end  of  execution  when  all 
data  in  the  dynamic  pool  is  written  into  the  appropriate  data  bases. 

Thus,  data  is  automatically  maintained  for  future  restarts. 

5.4  FINITE  System  Organization 
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The  logical  organization  of  FINITE  is  best  described  by  defining 
the  three  primary  capabilities  it  provides  structural  analysts.  These  are: 


1)  Definition  and  maintenance  of  tables  within  a data  base  that 
describe  the  characteristics  of  element  and  material  model 
library  routines; 

2)  Definition  and  maintenance  within  a data  base  of  geometric, 
topologic,  load  and  constraint  data  for  hierarchies  of  struc- 
tures ; 

3}  Computational  and  output  procedures  for  analyzing  structures 
and  reporting  the  results. 

To  support  the  system  capabilities  outlined  above  requires  a number 
of  processing  modules,  data  bases,  and  a mechanism  to  interface  between 
the  modules.  The  following  sections  briefly  describe  each  of  these 
conqjonents . 

5.4.1  Organization  of  Processing  Modules 

At  the  most  abstract  level,  the  physical  organization  of  FINITE 
consists  of  the  following  three  components: 

1)  The  POLO  supervisory  system; 

2)  FINITE  processing  modules  that  operate  under  control  of 
POLO  (POLO  subsystems); 

3)  Element  and  material  model  subprograms  initiated  by  FINITE 
subsystems . 

An  expanded  diagram  illustrating  additional  details  of  the  system 
organization  is  shown  in  Fig.  5.4.1.  POLO  is  represented  as  one  functional 
module  at  the  highest  level  with  the  primary  FINITE  subsystems  immediately 
below. 


Each  FINITE  processing  module  (currently  19  of  them)  is  a separate 
POLO  subsystem.  The  highest  level  processor,  FINITE,  ensures  existence 
of  appropriate  data  bases  and  translates  POL  commands  to  determine  which 
of  the  three  secondary  subsystems  to  initiate  (LIBRARY,  STORE,  or  COMPUTE). 
Each  of  these  subsystems  corresponds  directly  to  one  aspect  of  the  logical 
organization.  The  LIBRARY  subsystem  maintains  descriptions  of  elements 
and  material  models  supplied  by  their  developers.  STORE  translates  all 
POL  statements  defining  structural  geometries,  loads,  constraints,  etc. 
COMPUTE  contains  the  logic  required  to  perform  linear  analysis  by  initiating 
lower  level  modules  in  the  proper  sequence.  However,  COMPUTE  initiates 
NONCOMPUTE  to  direct  solution  of  nonlinear  structures.  As  shown  in  the 
figure,  some  modules,  for  example,  the  one  for  matrix  decomposition,  are 
invoked  by  several  modules  at  different  levels  in  the  hierarchy.  This 
illustrates  the  interdependency  of  processing  modules  but  also  shows 
that  they  are  self-contained  and  can  be  easily  initiated.  Lower  level 
subsystems  invoke  element  and  material  model  routines  to  compute  stiff- 
nesses, strains,  stresses,  etc. 

5.4.2  Data  Bases' 

The  conceptual  data  structure  for  finite  element  analysis  shown  in 
Fig.  5.2.2  has  been  partitioned  into  three  logical  POLO  data  bases  each 
of  which  resides  on  a separate  physical  file.  The  fourth  file  shown  in 
Fig.  5.4.1  is  the  POLO  system  data  base  that  contains  object  grammars  and 
file  definitions. 

The  data  base  denoted  LIBRARY  in  the  figure  contains  tables  describ- 
ing all  finite  elements  and  material  models  available  in  the  system.  Data 
in  the  tables  are  provided  by  developers  of  new  elements  and  material 
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models  through  POL  conmands  In  the  FINITE  language.  Geometric  properties 
of  common  structural  shapes  can  be  defined  In  the  library  to  minimize 
Input  data  during  structural  definition. 

The  primary  data  base  Is  denoted  WORKFILE  In  the  figure.  The 
Internal  form  of  all  structural  descriptions  and  computed  results  reside 
In  this  data  base.  Extensive  hlerarchlal  data  structures  similar  to  that 
shown  In  Fig.  5.3.3  are  defined  to  contain  element  and  substructure 
geometry,  topology,  stiffness,  loads,  displacements,  etc. 

A separate  data  base,  denoted  SOLVER,  Is  defined  for  solving  equations. 
Triangulated  equilibrium  equations  and  results  of  condensations  are 
retained  In  the  SOLVER  data  base.  These  are  available  to  process  addi- 
tional loading  conditions  as  well  as  to  retrieve  displacements  Inside 
condensed  substructures.  The  use  of  a data  base  separate  from  the 
WORKFILE  permits  optimal  allocation  of  the  equations  on  secondary  storage 
to  minimize  I/O  operations  during  equation  solving. 

5.4.3  Interfacing  Between  Subsystems 

During  analysis,  various  FINITE  subsystems  are  Invoked  by  other 
subsystems  to  perform  parts  of  the  solution  process.  Because  of  FINITE's 
flexibility,  the  exact  order  of  subsystem  Invocation  Is  highly  problem 
dependent  and  is  determined  by  the  system  as  execution  proceeds.  Thus, 
a flexible  technique  for  interfacing  between  subsystems  was  needed. 

Conceptually,  the  process  of  invoking  a hierarchy  of  subsystems  is 
Identical  to  a series  of  subprogram  calls  In  FORTRAN.  The  basic  actions 
required  to  initiate  a new  subsystem  are: 

1)  Save  the  status  of  any  variables  In  the  currently  executing 
subsystem  to  permit  resumption  of  processing  on 
the  lower  level  subsystem; 
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2)  Make  available  to  the  initiated  subsystem  controlling  parameters 
passed  by  the  initiating  subsystem; 

3)  Maintain  a hierarchial  list  of  subsystems  executed  to  reach 
the  current  level  to  enable  re-tracing  the  order  of  execution. 

The  first  two  functions  above  are  accomplished  through  a combination 
of  global  COMMON  parameters  and  a stack  of  "request"  vectors  maintained 
as  part  of  the  logical  data  space.  The  data  structure  for  this  is  shown 
in  Fig.  5.4.2.  The  third  function  above  is  performed  directly  by  the  POLO 
supervisor. 

To  initiate  a lower  level  module,  the  currently  executing  subsystem 
builds  a request  vector  of  parameters  for  initializing  the  subsequent 
module  and  pushes  the  request  onto  a stack.  The  global  parameter  REQLEVEL, 
available  to  all  subsystems,  indicates  the  current  stack  level.  The 
active  subsystem  then  requests  POLO  to  initiate  the  lower  level  sub- 
system. Once  executing,  the  new  subsystem  takes  its  instructions  from 
the  top  of  the  stack. 

This  simple  mechanism  is  used  in  all  modern  computers  (hardware), 
and  provided  FINITE  developers  with  an  ideal  scheme  for  maximum 
"appearance  of  integration"  while  maintaining  in  reality  19  different, 
and  loosely  interconnected  subsystems. 

5.5  Material  Modeling  in  FINITE 

The  unique  feature  of  nonlinear  material  modeling  in  FINITE  is  the 
isolation  of  elements  and  material  models  from  the  remainder  of  the  system 
by  standard  interfaces.  The  procedure  for  defining  new  linear  elements 
is  described  in  Ref.  [34].  Only  minor  changes  were  required  to  incorporate 
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additional  data  describing  nonlinear  elements.  The  discussion  here 
presents  the  new  aspects  of  FINITE  to  support  a general  nonlinear  material 
modeling  capability. 

The  Interface  mechanism  adopted  Is  Identical  for  all  material  models 
regardless  of  their  complexity.  An  Important  advantage  of  this  approach 
Is  that  developers  of  new  material  models  need  be  concerned  with  the 
organization  of  only  a small  “subspace"  of  the  entire  FINITE  system.  Once 
an  Individual  Is  familiar  with  the  Interface  scheme.  Implementation  of 
subsequent  material  models  Is  straightforward. 

In  addition  to  Isolating  material  models,  FINITE  performs  several 
"systems"  tasks  that  have  Inhibited  model  development  efforts  In  previous 
nonlinear  finite  element  systems.  Translation  of  material  model  Input 
data,  display  of  results,  memory  allocation,  data  retrieval,  and  compati- 
bility checking  between  elements  and  models  are  all  handled  by  FIMITE's 
material  processing  subsystems. 

The  following  section  describes  how  one  uses  existing  material 
models  In  FINITE  to  analyze  nonlinear  structures.  The  procedures  for  ; 

implementing  new  material  modeling  capabilities  are  described  in  section 
5.5.2. 

J 

5.5.1  Specification  of  Materials 

1 

During  the  Input  phase  of  analysis,  nonlinear  "materials"  are  | 

created  by  the  user  through  POL  statements  In  the  FINITE  language.  | 

Finite  elements  are  marked  materially  nonlinear  by  FINITE  Input  trans-  | 

lators  whenever  a nonlinear  material  Is  mentioned  as  part  of  an  element  j 

description.  | 
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A "material"  is  defined  by  selecting  three  possible  components 
from  finite's  library.  These  components  are: 

1)  A material  model; 

2)  A stress-strain  function  compatible  with  the  material  model 
(optional); 

3)  An  initial-strain  function  compatible  with  the  material  model 
(optional). 

The  material  model  determines  stresses  in  the  material  given  the 
previous  strain  and  stress,  new  strain,  loading  history,  etc.  Most 
classical  models  based  on  plasticity  theory  utilize  a yield  or  loading 
surface  with  some  flow  rule.  However,  FINITE's  material  system  is  in 
no  way  restricted  to  plasticity  type  behavior.  The  stress-strain 
function  provides  the  selected  model  with  physical  properties  of  the 
material,  for  example,  results  of  a uniaxial  tension  test.  The  initial- 
strain  function  supplies  the  model  with  time  dependent  effects  on  material 
behavior. 

The  following  POL  statements  define  a simple  material: 

MATERIAL  A36-STEEL  TYPE  VON-MISES 
PROPERTIES  MAXINCR  30  ALPHA  0.05  DIVERGE  500 
USE  STRESS-STRAIN  FUNCTION  SEGMENTAL 
PROPERTIES  E 30000.  NU  0.3  SEGMENTAL, 

NUMPOINTS  3, 

STRAINS  0.0  0.001,  0.008, 

STRESSES  0.  36.  42. 

The  material  name  is  A36-STEEL  and  its  behavior  is  modeled  by  the 
VON-MISES  material  model  selected  from  FINITE's  library.  Users  may 
reference  A36-STEEL  when  defining  the  properties  of  finite  elements. 

PROPERTIES  are  parameters  defi..ed  for  a model  or  function  for  which 
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users  can  supply  values  during  analysis.  In  the  above,  ALPHA,  MAXINCR, 
and  DIVERGE  control  computational  algorithms  within  the  VON-MISES 
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material  model. 

The  stress-strain  function  SEGMENTAL  selected  from  the  library 
approximates  the  uniaxial  tension,  stress-strain  curve  by  a sequence  of 
straight  lines.  Data  describing  the  curve  are  also  given  through  a 
PROPERTIES  statement  following  the  function  declaration. 

The  input  data  shown  above  is  given  prior  to  requesting  a nonlinear 
analysis.  Once  the  incremental-iterative  solution  process  has  begun,  a 
user  can  modify  values  for  model  and  function  PROPERTIES  whenever  he  has 
control  between  load  steps.  The  existing  material  is  referenced  with  the 
word  ACCESS  given  prior  to  MATERIAL.  Only  those  model  or  function 
properties  that  have  changed  must  be  re-entered.  Different  functions 
may  also  be  associated  with  the  model.  This  feature  is  most  convenient 
for  altering  tolerances,  model  algorithms,  etc.  as  the  nonlinear  solution 
progresses. 

5.5.2  Entry  of  New  Models  and  Functions 

The  definition  of  new  material  models  and  functions  is  divided 
into  two  parts  — tables  and  corresponding  subroutines.  Information 
stored  in  the  tables  is  used  by  FINITE  to  make  compatibility  checks, 
translate  input  data  and  allocate  memory  prior  to  initiating  the  subroutines. 
The  individual  defining  a new  model  or  function  enters  both  tables  and 
subroutines.  The  table  definition  process  is  termed  "soft"  because  the 
developer  may  alter  the  information  at  any  time  by  simply  re-entering 
the  table  data  during  a FINITE  run.  This  is  convenient  for  developers 
working  with  experimental  material  models  that  change  frequently. 

The  POL  statements  defining  the  VON  MISES  table  data  are  listed  below: 


DEFINE  MATERIAL  MODEL  VON-MISES 
SUBROUTINE  1 

MATERIAL  STRESS-STRAIN  FUNCTIONS 

SEGMENTAL,  RAMBERG-OSGOOD , ELASTO- PLASTIC 
HISTORY  5 WORDS 
ELEMENTS 

CSTRIANGLE,  L2DISOP,  Q2DISOP 

PROPERTIES 

PSTRAIN  LOGICAL  FALSE 
TOLERANCE  REAL  O.OOl 


END  OF  PROPS 
SYMMETRY 
END  OF  MODEL 

The  SUBROUTINE  statement  associates  the  tables  generated  via  the  above 
commands  with  FORTRAN  subroutines  that  perform  material  model  computations. 

Names  of  functions  and  elements  compatible  with  the  material  model 
supplied  in  the  tables  permit  the  FINITE  input  translators  to  check 
compatibility  of  the  material  model,  stress-strain  and  initial-strain 
functions,  and  the  element  type  during  the  input  phase  of  analysis.  By 
delaying  all  consistency  checks  until  an  analysis  is  performed,  developers 
may  define  elements,  models,  and  functions  in  the  library  in  any  desired 
order,  i.e.,  reference  to  non-existent  elements  and  functions  is  per- 
mitted. 

The  HISTORY  statement  declares  the  amount  of  working  storage  space 
to  reserve  for  each  strain  point  of  each  element  that  references  the 
material  model.  Information  stored  in  this  space  is  made  available  to 
the  model  each  time  the  strain  point  is  processed  during  nonlinear  analysis. 
Material  models  may  use  the  space  to  retain  load  path  dependent  parameters 
governing  material  behavior  etc.  ...  FINITE  is  unaware  of  the  type  or  mode 
of  data  within  the  history  space. 


PROPERTIES  permit  material  model  developers  to  define  keywords  and 
default  values  for  parameters  associated  with  the  model.  Users  may 
optionally  override  the  default  value  declared  in  the  table  during 
problem  solution  as  shown  previously.  Real,  integer,  logical,  and 
alphanumeric  modes  of  data  are  acceptable  as  properties.  Properties  may 
be  either  scalars  (single-valued)  or  vectors  (multi-valued). 

Most  material  models  generate  a symmetric  matrix  relating  incremental 
changes  of  stress  and  strain.  The  SYMMETRY  statement  declares  this 
result.  Models  that  compute  a nonsymmetrical  matrix,  for  example,  those 
using  non-associated  flow  rules,  must  declare  this  with  a NONSYMMETRY 
statement.  FINITE  automatically  marks  elements  associated  with  non- 
symmetric  materials  as  having  nonsymmetric  stiffness  matrices.  Similarly, 
substructures  in  which  these  elements  appear  are  also  nonsymmetric  as  are 
the  final  equilibrium  equations.  FINITE  handles  this  situation  without 
user  or  developer  intervention  from  input  translation  through  solution 
of  the  nonsymmetric  equations. 

Tables  defining  stress-strain  and  initial-strain  functions  are 
simpler  than  those  for  a material  model.  The  SEGMENTAL  function  referred 
to  in  the  previous  example  is  defined  by  the  following: 

DEFINE  MATERIAL  STRESS-STRAIN  FUNCTION  SEGMENTAL 
SUBROUTINE  1 
PROPERTIES 
YMODULUS  REAL  0.0 
NUMPOINTS  INTEGER  0 
STRAINS  VECTOR  REAL  0.0 
STRESSES  VECTOR  REAL  0.0 


ENO  OF  PROPERTIES 
END  OF  FUNCTION 
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Complete  table  Input  data  for  several  material  models  and  functions  are 
given  in  Appendix  B. 

5.5.3  Model  and  Function  Subprograms  ^ 

During  problem  solution,  the  FINITE  material  processor  calls  material 
model  subprograms  to  perform  computations.  Based  on  the  type  of  calcula- 
tions required,  the  processor  locates  and  makes  present  in  memory  all 
data  needed  by  material  models  and  functions.  Data  for  a single  strain 
point  is  extracted  and  passed  to  the  material  model  for  computation. 

This  procedure  continues  until  all  strain  points  of  materially  nonlinear 
elements  of  the  structure  (and  substructures)  have  been  processed. 

Names  of  material  model  routines  are  MTMXX  where  XX  is  the  two-digit 
sequence  containing  the  SUBROUTINE  number  and  is  taken  from  the  library 
tables.  The  calling  sequence  is 

SUBROUTINE  MTMXX  (PROPPT, PROPS, HISTRY,DMTRIX,NSTRM,OLDEPS, 
OLIEPS  ,DIEPS ,0LDSI6 ,NEWSI6 ,NWIEPS ,6ENDAT , 

TRANS) 

A summary  of  the  most  important  parameters  is  given  below.  Those  narked 
with  an  (*)  are  updated  by  the  model  depending  on  the  type  of  computation 
requested  by  FINITE. 


PROPS  - is  the  vector  of  model  property  values. 

HISTRY  - (*)  is  the  history  data  vector  for  the  strain  point. 

DMTRIX  - (*)  is  the  tangent  modulus  matrix  for  the  strain  point 
currently  incorporated  in  the  element  and  structure  stiffness 
matrix. 

OLDEPS  - is  the  vector  of  total  strains  for  the  previous  load 
increment  or  iteration. 
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OLIEPS  - is  the  vector  of  initial  strains  included  in  OLOEPS,  i.e., 
OLDEPS-OLIEPS  = total  mechanical  strain. 

DEPS  - is  the  vector  of  incremental  strains  for  this  load  increment 
or  iteration. 

DIEPS  - (*)  is  the  vector  of  initial  strains  included  in  the  incre- 
mental strain  increment  (DEPS).  For  a load  increment  (iteration 
* 1),  this  vector  contains  initial  strains  due  to  the  real 
applied  incremental  load.  Thereafter,  it  contains  whatever 
was  entered  in  the  NWIEPS  vector  during  the  previous  call  to 
the  material  model,  i.e.,  incremental-iterative  creep,  shrinkage, 
etc.,  strains. 


OLDSIG  - is  the  vector  of  total  stresses  for  the  previous  load 
increment  or  iteration. 


NWIEPS  - (*)  is  a zeroed  vector  in  which  the  model  may  place  new 
incremental  initial  strains  due  to  creep,  shrinkage,  etc.,  that 
occur  during  a load  step  or  iteration.  These  initial  strains 
are  converted  into  equivalent  forces  by  the  element  residual 
load  routines. 

GENDAT  - (*)  is  a vector  of  general  data  some  of  which  is  updated 
by  the  model  to  indicate  results  of  computation. 

TRANS  - a vector  of  geometric  parameters  generated  by  the  element 
strain  routine  for  use  in  transforming  from  one  definition  of 
strain-stress  to  another  for  geometrically  nonlinear  problems. 
Generally  this  vector  contains  data  describing  geometry  of 
deformed  elements. 

Material  model  subprograms  reference  stress-strain  or  initial -strain 
functions  via  standard  FINITE  system  interfaces  termed  MTMDSS  and  MTMDIS. 
Material  models  communicate  with  functions  through  a set  of  parameters 
termed  a communication  vector.  The  protocol  for  data  in  the  vector  is 
established  by  the  particular  function  used.  Most  functions  require  and 
return  similar  data,  thus  vector  protocols  are  nearly  identical.  To 
reference  a function,  the  material  model  builds  the  communication  vector 
and  issues  a call  to  MTMDSS  or  MTMDIS.  These  FINITE  routines  invoke  the 
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communications  vector  after  which  successive  returns  are  made  to  the  material 
model . 

It  should  be  restated  that  the  use  of  functions  is  not  mandatory. 
Material  models  that  require  very  simple  material  characteristics  may 
include  them  in  the  model  properties  and  thus  eliminate  the  need  for 
functions. 

Names  of  function  subroutines  are  MTSXX  and  MTIXX  for  stress-strain 
and  initial-strain  functions,  respectively.  XX  is  the  subroutine  number 
declared  for  the  function  in  the  library.  The  calling  sequences  are 

SUBROUTINE  MT IXX ( COMVEC , PROPPT , PROPS .GENDAT ) 
and 

SUBROUTINE  MTSXX (COMVEC .PROPPT .PROPS .GENDAT) 

The  parameters  are: 

COMVEC  - the  cortmuni cation  vector  constructed  by  the 
material  model  routine; 

PROPPT  - a vector  of  subscripts  used  to  locate  property 
values  in  PROPS; 

PROPS  - the  vector  of  function  property  values; 

GENDAT  - a vector  of  general  information  passed  to  the 
function.  / 

/ 

5.6  Nonlinear  Processors  of  FINITE 

After  the  input  phase  of  analysis,  users  may  request  computation 
of  nonlinear  results  for  any  number  of  consecutive  load  steps.  Sub- 
system NONCOMPUTE  contains  the  logic  to  perform  the  nonlinear  solution 
given  the  I'ser's  request  in  a suitable  internal  form. 

Fig.  5.6.1  is  a simplified  block  diagram  showing  the  major  flow  of 
control  through  NONCOMPUTE.  A special  data  base  initialization  is 
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required  If  the  request  Is  for  computation  of  load  step  number  one. 

Another  separate  section  of  the  NONCOMPUTE  subsystem  performs  the  detailed 
sequence  of  operations  necessary  to  effect  solution  for  a single  load 
step. 

After  solution  of  all  requested  load  steps,  control  Is  returned  to 
the  FINITE  subsystem.  If  tn  an  Iteractive  environment,  the  user  may  ask 
for  output,  modify  solution  parameters  and  continue  the  analysis,  etc. 

For  batch  operation,  the  next  corrmands  In  the  Input  stream  are  processed. 
Eventually,  data  bases  are  secured  In  preparation  for  future  analysis 
restart  and  execution  Is  terminated. 

The  two  major  components  of  NONCOMPUTE,  data  structure  Initialization 
and  details  of  a single  step  solution,  are  discussed  In  the  following 
sections. 

5.6.1  Data  Structure  Initialization 

A convenient  feature  of  FINITE  Is  that  nonlinear  structural  hierarchies 
are  defined  by  users  exactly  as  for  linear  analysis.  In  the  linear 
version  of  FINITE,  data  structures  naturally  suited  to  take  advantage  of 
Identical  elements  and  substructures  were  developed.  The  stiffness 
matrix  and  equivalent  nodal  loads  for  a substructure  are  computed  just 
once  and  are  referred  to  by  each  occurrence  of  the  substructure  when  used 
as  an  element  In  a higher  level  structure.  This  procedure  Is  not  valid 
for  nonlinear  substructures  because  they  have  unique  stiffness  matrices 
and  loads  for  each  occurrence  In  the  hierarchy.  An  "expansion"  process 
was  designed  and  Implemented  that  makes  each  occurrence  of  a nonlinear 
element  or  substructure  a unique  component  In  the  structural  hierarchy. 
Expansion  Is  totally  transparent  to  the  user  and  does  not  destroy  the 
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original  linear  data  structure  thereby  permitting  a linear  analysis  of 
the  structure  if  desired  at  a later  time. 

Expansion  of  the  structural  hierarchy  is  performed  during  solution 
for  the  first  load  step.  Important  components  of  the  process  are  out- 
lined in  Fig.  5.6.2.  The  order  of  operations  is  designed  to  eliminate 
wasteful  computation  of  stiffness  matrices  for  multiple  occurrences  of 
nonlinear  elements  and  substructures  which  are  always  identical  (linear) 
in  the  first  load  step.  Once  the  linear  stiffness,  based  on  the  unexpanded 
hierarchy,  has  been  computed,  nonlinear  elements  and  substructures  are 
duplicated  in  the  hierarchy  for  each  occurrence.  Data  independent  of 
nonlinear  effects,  for  example  incidences  and  coordinates,  are  not 
duplicated;  rather  data  for  the  original  occurrence  is  referenced.  After 
expansion  of  the  hierarchy,  additional  data  structures  are  initialized 
for  nonlinear  data.  Space  is  created  for  the  incremental  constitutive 
matrices  and  history  parameters  for  all  materially  nonlinear  elements. 

For  geometrically  nonlinear  elements,  a data  structure  is  initialized  to 
contain  element  local  displacements  for  rapid  access  independent  of  the 
structure  in  which  the  elements  reside. 

5.6.2  Solution  for  a Load  Step 

The  order  of  computations  for  a single  load  step  is  shown  in  Fig. 
5.6.3.  At  the  NONCOMPUTE  level,  the  solution  process  is  deceptively 
simple.  The  load  step  driver  of  NONCOMPUTE  invokes  appropriate  lower 
level  subsystems  in  the  proper  sequence  corresponding  to  the  solution 
parameters  declared  by  the  user.  Complexities  due  to  multi-level  sub- 
structuring and  condensation  are  handled  by  the  individual  lower  level 
subsystems . 
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The  first  Iteration  of  a load  step  corresponds  to  application  of 
the  real  load  Increment.  Subsequent  Iterations,  If  permitted  and  required, 
are  performed  at  a constant  real  load  level  to  remove  the  residual  nodal 
loads  and  to  bring  the  structure  Into  an  equilibrium  configuration. 

An  Important  aspect  of  the  solution  process  as  outlined  In  the  flugre 
Is  data  security  in  the  event  of  analysis  failure.  If  the  solution  for  a 
load  step  falls  to  converge,  the  user  can  redefine  the  load  step  and 
request  a new  solution.  To  permit  this,  the  system  retains  enough  data 
between  load  steps  to  reconstruct  the  solution  status  at  any  point. 

Limits  on  the  amount  of  secondary  storage  available  prevent  saving  of 
previous  tangent  stiffness  matrices  and  the  triangulated  equilibrium 
equations  for  each  step.  Only  data  necessary  to  reconstruct  element  and 
structure  stiffnesses  are  saved.  For  example.  If  material  nonl Inearl  ties 
are  present,  the  Incremental  constitutive  matrices  and  history  parameters 
are  retained.  Thus,  If  a stiffness  regeneration  Is  required  before 
processing  can  resume,  the  user  simply  Indicates  this  through  a solution 
parameter  then  requests  solution  for  the  same  step  again.  Nonlinear 
FINITE  subsystems  automatically  destroy  and  recompute  data  rendered 
Invalid  as  a result  of  the  re-analysis. 

For  large  structures  or  structures  analyzed  with  many  load  steps, 
the  amount  of  data  on  secondary  storage  may  become  very  large.  Commands 
are  available  for  users  to  Indicate  that  results  for  specific  load  steps 
need  no  longer  be  retained.  Space  on  secondary  storage  occupied  by  this 
data  can  then  be  reused. 
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In  the  current  nonlinear  version,  once  computations  begin  for  a 
load  step  the  user  cannot  regain  control  over  execution  until  the  step 
computations  are  completed.  In  the  future,  provisions  should  be  made  to 
permit  resumption  of  load  step  computations  in  the  event  of  nonconvergence. 


CHAPTER  6 


NUrCRICAL  EXAMPLES 

6. 1 General 

The  nonlinear  capabilities  in  the  FINITE  system,  as  designed  and 
Implemented  In  this  study,  are  demonstrated  through  a series  of  example 
analyses  In  this  chapter.  A variety  of  structural  configurations  and 
nonlinear  effects  are  considered.  Realistic  structures  In  both  size  and 
complexity  were  chosen  for  examples  to  Illustrate  the  practical  aspects 
of  the  system. 

These  problems  are  not  intended  to  examine  the  adequacy  of  any 
particular  finite  element  or  material  model  formulations.  They  demon- 
strate that  a)  material  models  can  be  easily  Implemented,  b)  FINITE  can 
handle  a variety  of  nonlinear  behaviors,  and  c)  and  the  system  Is  relatively 
easy  to  use.  Elements  and  models  employed  In  the  examples  have  already 
proven  effective  In  other  studies.  Emphasis  here  Is  placed  on  structural 
modeling,  data  generation,  ease  of  input,  multiple  analysis  restarts  and 
other  aspects  of  the  system  experienced  directly  by  users.  Of  particular 
Interest  Is  the  application  of  substructuring  and  condensation  to  reduce 
analysis  cost. 

The  first  two  examples  Illustrate  the  analysis  of  structures  composed 
of  long  slender  members  In  which  nonlinear  behavior  occurs  as  a result  of 
large  rigid  body  rotations.  A space  truss  dome  and  plane  frame  building 
are  the  two  structures  considered.  In  both  problems  the  Incremental 
solution  method  provides  deformations  prior  to  buckling  as  well  as  the 
approximate  critical  load  for  the  structure. 


The  last  three  examples  consider  the  effects  of  nonlinear  material 
behavior  in  solid  mechanics  problems.  Plastic  deformation  near  a junction 
in  a steel  axisymmetric  pressure  vessel  using  the  Von  Mises  criterion  is 
examined  in  the  first  solid  mechanics  problem.  The  second  problem 
simulates  a tunnel  opening  in  rock  using  the  Orucker-Prager  yield  criterion. 
The  last  example  is  a three  dimensional  plastic  analysis  of  a perforated 
thick  plate  subjected  to  uniform  pressure. 

6.2  Space  Truss 

Reticulated  domes  constructed  of  truss- type  members  have  long  been 
an  economical  solution  for  roof  structures.  Most  analyses  in  the  past 
were  based  on  linear  theory  even  though  shallow  domes  were  known  to  respond 
nonlinearly  under  design  loads  due  to  changes  in  geometry.  Under  higher 
load  levels,  shallow  domes  exhibit  snap-through  buckling  behavior.  Even 
more  important  from  a design  consideration  is  the  sensitivity  of  these 
structures  to  nonsymmetric  loads  resulting  from  wind  pressure  and  snow. 

One  particularly  interesting  geometric  configuration  is  the  Schwedler 
dome  shown  in  Fig,  6,2.1.  There  are  24  equal  divisions  on  each  latitudinal 
circle  and  4 equal  divisions  along  each  meridian.  All  nodes  of  the  dome 
lie  on  the  surface  of  a sphere.  Base  nodes  are  completely  restrained 
against  translation.  Dimensions  and  material  properties  are  shown  on  the 
figure. 

The  structure  has  cyclic  symmetry  but  does  not  have  radial  symmetry; 
thus  all  nodes  must  be  included  in  any  analysis.  Under  a uniform  vertical 
load  the  structure  twists  counter-clockwise  about  the  vertical  axis. 

The  finite  element  formulation  for  a general  axial  force  member 
including  all  nonlinear  strain  terms  is  given  in  Appendix  A.  The  significance 
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of  retaining  all  nonlinear  terms  has  been  demonstrated  In  Ref.  [30] 
which  also  solved  the  Schwedler  dome  structure.  Omission  of  certain  terms 
leads  to  straining  of  an  element  subjected  to  a rigid  body  rotation. 

The  analysis  here  considers  two  separate  loading  conditions.  First, 
a uniform  gravity  load  of  60  psf  Is  applied  over  the  entire  dome.  The 
second  consists  of  a uniform  vertical  load  acting  on  one  half  of  the  dome 
designated  ABC  In  the  figure.  The  magnitude  of  this  load  is  Increased 
In  a series  of  Increments  until  Instability  occurs.  For  both  loading 
conditions,  equivalent  nodal  loads  corresponding  to  an  applied  vertical 
load  of  1 psf  form  the  basis  of  nonlinear  load  Increments. 

Two  structural  models  were  developed  for  analysis.  The  first  model 
Includes  all  nodes  and  elements  as  shown  in  Fig.  6.2.1  and  is  termed  the 
standard  model.  An  examination  of  the  structure  Indicates  that  nonlinear 
effects  should  be  significant  In  the  nearly  flat  part  of  the  dome  shown 
In  Fig.  6.2.2c.  Steeper  meridians  In  the  lower  portion  are  not  signifi- 
cantly Influenced  by  small  geometry  changes.  The  lower  portion  Is  modeled 
as  a linear  substructure  with  condensation  applied  to  eliminate  all  nodes 
except  those  connecting  the  top  ring.  The  substructured  model  Is  shown 
In  Fig.  6.2.2, 

The  Incremental -Iterative  analysis  procedure  was  carried  out  using 
two  load  steps  for  the  full  gravity  loading.  Five  load  steps  were  used 
to  attain  the  Instability  load  for  the  partial  loading.  Iterative 
corrections  for  residual  loads  were  performed  at  each  step.  Convergence 
checks  were  based  on  the  residual  load  vector  norm.  A rapidly  Increasing 
residual  load  norm  Indicates  the  Instability  load  has  been  reached.  At 
the  critical  load  level,  the  structure  attempts  to  snap-through  to  the 
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alternate  equilibrium  configuration  Since  this  configuration  is  of  no 
practical  interest,  the  analysis  was  terminated. 

Portions  of  the  FINITE  input  data  for  the  two  structural  models 
are  listed  in  Fig.  6.2.3  and  6.2.4.  Both  sets  of  input  data  illustrate 
the  problem-oriented-language  commands  to  describe  a structure,  its  loads, 
and  the  analysis  procedure.  Input  data  for  both  models  is  straight- 
forward and  self- documenting. 

Results  for  the  uniform  load  case  are  shown  in  Fig.  6.2.5.  The  numer- 
ical values  given  are  for  the  standard  model.  Deflections  for  the  substruc 
tured  model  (not  shown)  are  slightly  larger  as  the  stiffening  effect  due  to 
geometry  changes  is  not  fully  accounted  for.  However,  the  differences 
are  less  than  3%.  Results  shown  in  Fig.  6.2.6  for  the  partial  loading 
condition  are  for  the  load  level  (29  psf)  just  prior  to  instability.  The 
standard  model  did  not  converge  at  a load  level  of  29.5  psf.  The  sub- 
structured model  did  not  converge  at  a load  level  of  30  psf.  Little 
difference  is  noted  between  the  standard  and  substructured  models  thus 
verifying  the  original  assumption  regarding  the  location  of  nonlinear 
effects.  The  primary  difference  between  the  standard  and  substructured 
models  is  the  computer  time  required  for  analysis.  For  the  full  uniform 
loading  condition,  the  substructured  model  required  745i  as  much  time  as 
the  standard  model  (2  load  steps).  More  noticeable  savings  are  achieved 
for  the  partial  loading  condition.  For  that  case,  the  substructured  model 
used  only  53%  as  much  time  as  the  standard  model  (5  load  steps).  Even 
greater  savings  would  occur  with  a more  detailed  analysis  involving  a 
larger  number  of  steps. 
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6.3  Second  Order  Frame  Analysis 

Linear  analysis  of  plane  frame  structures  neglects  the  interaction 
of  axial  forces  and  bending  stiffness.  This  interaction  becomes  signifi- 
cant in  braced  frames  with  large  column  loads  and,  in  all  unbraced 
frames.  Current  design  codes  require  amplivi cation  of  member  forces 
computed  by  a linear  analysis  to  account  for  interaction  effects. 
Magnification  factors  applied  to  linear  results  have  been  developed  from 
approximate  nonlinear  analyses.  For  practical  design,  stresses  remain 
well  below  the  proportional  limit;  thus,  the  determination  of  axial  load 
effects  on  bending  action  involves  only  geometrical  nonlinearity. 

Most  design  codes  permit  a more  accurate  determination  of  member 
forces  in  frame  structures  by  a "second  order"  analysis  procedure.  A 
second  order  analysis  must  rationally  determine  the  effect  of  axial  loads 
on  rotational  stiffness  coefficients  in  braced  frames  and  the  P-delta 
effect  in  unbraced  frames. 

Large  compressive  column  loads  in  braced  frames  reduce  the  moment 
exerted  by  the  column  on  any  rigidly  framed  beams,  thereby  increasing 
deflections  and  positive  beam  moments.  Buckling  of  a braced  frame  occurs 
when  the  applied  load  drastically  reduces  column  rotational  stiffness 
coefficients. 

Vertical  loads  acting  through  sway  deflections  of  unbraced  frames 
produce  additional  bending  moments  and  shears.  This  phenomena  is  termed 
the  P-delta  effect.  Member  forces  in  unbraced  frames  determined  by  a 
second  order  analysis  include  this  effect  and  are  not  modified  by  the 
amplification  factors  in  design  codes.  Unbraced  frame  members  are 
proportioned  to  resist  second  order  forces  using  the  same  procedures  for 
braced  frames. 
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I The  nonlinear  plane  frame  element  in  FINITE  has  three  degrees  of 

\ freedom  (U,  V,  0^)  at  each  node  as  shown  in  Fig.  6.3.1.  Element  stiffness 

I characteristics  are  determined  from  the  classical  differential  equation 

' of  the  elastic  curve  including  axial  load 


E I V"(X)  = (±)  P V(X)  (6.3.1) 

where  E is  Young's  modulus,  I the  moment  of  inertia  in  the  plane  of 
bending,  P the  axial  force,  and  V(x)  is  the  lateral  deflection  along  the 
member  measured  from  the  undeformed  configuration.  Tensile  and  compres- 
sive axial  forces  are  permitted.  The  element  secant  stiffness  is  derived 
in  Refs.  [9,  39]  using  stability  functions  obtained  from  homogeneous 
solutions  of  this  equation. 

Approximations  employed  in  the  development  of  the  differential  equation 
must  be  carefully  delineated  for  practical  nonlinear  analysis.  They  are: 

1)  The  curvature  of  the  member  is  infinitesimal  and  may  be  ade- 
quately described  by  the  second  derivative  of  lateral  displace- 
ment, V"(X).  This  follows  directly  by  assuming  the  member 
centerline  slope  squared,  [V'(X)], is  negligible  compared  to  1. 

2)  Axial  forces  in  members  may  be  determined  by  the  equation 

P = E A e/L  (6.3.2) 


where  e is  the  member  elongation  given  by  as  shown 

in  Fig.  6.3.1.  A is  the  cross-sectional  area  and  L is  the 
undeformed  length.  This  assumption  implies  infinitesimal 
member  length  changes  and  that  changes  in  direction  of  member 
end  forces  due  to  the  sway  angle  are  neglected. 
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3)  The  deflected  shape  of  a laterally  loaded  menter  subjected 


to  an  axial  force  Is  adequately  described  by  equivalent  forces 
applied  at  the  ends  as  determined  from  the  deflected  shape 
without  axial  force. 

These  assumptions  have  no  significant  consequence  In  practical 
analysis  and  design.  Drift  limitations  Imposed  by  building  codes  render 
1)  and  2)  valid.  Normally  beams  are  the  only  members  subjected  to  lateral 
loads  and  the  effect  of  small  axial  forces  Is  very  insignificant;  thus, 
the  third  approximation  Is  valid. 

The  nonlinear  plane  frame  element  In  FINITE  adequately  predicts  the 
second  order  effects  as  required  by  design  codes  for  use  as  an  alternate 
analysis  method.  The  element  predicts  classical  linear  buckling  loads 
for  braced  frames  and  for  unbraced  frames  If  the  pre-buckling  deformations 
are  small.  The  following  two  examples  solved  with  FINITE  Illustrate 
simple,  but  practical  applications  of  the  nonlinear  capability. 

P-Delta  Analysis 

An  unbraced  plane  frame  structure  with  dimensions  and  loads  shown 
in  Fig.  6.3.2  was  analyzed  for  P-delta  effects.  Deflections  from  the 
ficticious  lateral  load  method  described  by  Adams  [1,2]  are  compared 
with  FINITE  results.  Following  Adams,  the  column  loads  for  P-delta 
analysis  are  Increased  by  1.7  over  working  loads  to  approximate  ultimate 
conditions.  Although  some  researchers  disagree  with  this  philosophy, 
the  Increased  loads  are  used  here  to  permit  direct  comparison  with  the 
published  results  of  Adams.  Fig.  6.3.3  summarizes  the  lateral  displace- 
ments computed  by  a standard  linear  analysis,  Adam's  method,  and  FINITE. 
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Displacements  computed  by  FINITE  are  larger  than  those  predicted  by  the 
lateral  force  method.  This  follows  because  the  shape  functions  used  in  the 
FINITE  formulation  are  exact  within  classical  theory.  The  lateral  force 
method  tries  to  correct  for  unbalanced  element  end  forces  using  the 
linear  stiffness  matrix.  Since  both  models  are  "compatible**  in  finite 
element  terminology,  the  one  based  on  correct  shape  functions  is  less 
stiff. 

Input  data  required  to  analyze  the  structure  with  FINITE  is  listed 
in  Fig.  6.3.4.  All  beams  were  assumed  to  respond  linearly.  Two  load 
steps  were  used;  each  step  required  3 iterations  for  convergence  with 
stiffness  updates  performed  before  the  second  iteration  of  each  step. 


Linear  Buckling  Analysis 

The  nonlinear  plane  frame  element  in  FINITE  can  be  used  to  compute 
approximate  buckling  loads  as  illustrated  in  this  example.  The  plane 
frame  structure  shown  in  Fig.  6.3.5  is  subjected  to  a distribution  of 
vertical  loads  corresponding  to  self  weight,  superimposed  dead  load,  and 
uniform  live  load  on  all  girders.  Very  small  lateral  loads  insure  the 
frame  will  buckle  in  a sidesway  mode.  The  initial  loading  pattern  repre- 
sents the  magnitudes  for  a unit  value  of  the  loading  parameter  x.  During 
analysis,  proportional  loads  are  applied  by  simply  incrementing  X. 
Instability  of  the  frame  occurs  when  lateral  displacements  increase 
rapidly  for  small  changes  in  x.  The  critical  load  level  is  denoted  x^^. 

Estimates  of  x^^  based  on  displacements  for  value  of  ^ 
computed  from 
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where  and  are  the  linear  and  nonlinear  displacements  of  any  story  ; 

corresponding  to  loading  level  x.  A derivation  and  discussion  of  this 

! ; 

useful  formula  are  given  In  Ref.  [39].  An  average  x^^  computed  from  x^^ 
of  each  story  provides  a guideline  for  the  next  Increment  of  loading  to 
apply. 

Fig.  6.3.6  presents  a summary  of  the  analysis  results.  The 
P-delta  load  deflection  curve  for  the  top  story  Is  given  In  Fig.  6.3.7. 

i 

6.4  Axisymmetric  Pressure  Vessel 

A steel  axisymmetric  pressure  vessel  with  a cylinder  to  sphere 
junction  Is  shown  In  Fig.  6.4.1.  Design  of  these  vessels  normally  follows 
elastic  theory  supplemented  by  standard  design  codes  of  practice.  Stress 
distributions  are  often  determined  from  simple  membrane  theory  with  stress 
concentration  factors  to  account  for  bending  and  localized  effects.  The 
concentration  factors  are  generally  the  result  of  elastic  thin  shell 
solutions;  however,  shell  theories  are  unable  to  easily  account  for  the 
stiffening  effect  of  the  junction  weld.  Upper  and  lower  bounds  for  the  | 

critical  pressure  can  be  determined  by  approximate  limit  analyses. 

i 

The  finite  element  method  is  naturally  suited  for  analysis  of  this  1 

type  structure.  Nonlinear  behavior  Is  Isolated  within  a small  distance  | 

i 

of  the  junction.  This  example  Illustrates  1)  the  specification  of 
materially  nonlinear  problems  In  FINITE  and  2)  substructuring  and  conden- 
sation In  stress  concentration  problems. 

A vessel  tested  by  DInno-GIll  [16]  was  chosen  for  analysis  to  permit 

I 

comparison  of  the  finite  element  solution  with  experimental  results.  i 


Geometry  of  the  vessel  Is  shown  in  Fig.  6.4.1.  Failure  occurs  well  below 
the  level  of  loading  to  cause  large  deformations;  thus  only  material 
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nonlinearity  is  considered.  The  vessel  was  constructed  of  mild  steel 
and  stress  relieved  after  completion  of  all  welding.  The  Von  Mises 
yield  criterion  adequately  predicts  material  behavior  under  these  con- 
ditions. 

Eight  node  isoparametric  elements  were  chosen  to  idealize  the 
vessel.  Strains  vary  linearly  over  each  element.  Nonlinear  elements  are 
used  only  in  the  immediate  vicinity  of  the  cylinder-sphere  junction. 

Remaining  elements  in  the  sphere  and  cylinder  are  assumed  to  respond 
linearly  for  all  load  levels.  The  finite  element  mesh  is  shown  in  Fig. 

6.4.2.  Reduced  integration,  2x2  rule,  employed  in  the  linear  regions  ^ 

reduces  cost  and  improves  element  performance.  In  the  nonlinear  elements,  I 

full  3x3  integration  permits  early  detection  of  yielding  due  to  high  ] 

bending  stresses  near  the  inner  and  outer  vessel  wall  surfaces. 

The  multi-level  substructuring  features  of  FINITE  employed  in  this 
problem  reduce  input  data  and  analysis  costs.  The  linear  portion  of  the 
cylinder  is  modeled  with  two  levels  of  substructuring  to  take  advantage 
of  identical  elements  in  the  cylinder  wall.  Static  condensation  applied  J 

to  both  the  linear  sphere  and  cylinder  substructures  eliminates  all  but 
three  nodes  necessary  to  connect  with  nonlinear  elements  as  shown  in  the 

I 

figure.  The  final  structure  consists  of  the  36  nonlinear  elements  and 
the  two  condensed  substructures.  A unit  pressure  load  applied  to  the 
inside  face  of  all  elements  in  the  linear  substructures  and  nonlinear 
elements  constitutes  the  loading  condition.  Scalar  multiples  of  this 
loading  define  the  individual  step  load  increments.  Other  aspects  of  the 
analysis  are  described  in  the  input  data  shown  in  Fig.  6.4.3. 
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The  spread  of  plastic  zones  for  the  levels  of  Internal  pressure 
are  shown  In  Fig.  6.4.4.  The  load  deflection  curve  for  node  59,  shown 
In  Fig.  6.4.5,  agrees  well  with  experimental  results. 

Several  aspects  of  the  solution  process  from  the  user  viewpoint  are 
worthwhile  considering.  Data  generation  comnands  simplify  the  Input  of 
nodal  coordinates  and  element  topology.  Loadings  on  lower  level  sub- 
structures are  carried  upward  through  the  hierarchy  via  the  EXTERNAL 
LOADS  command.  The  nonlinear  processor  of  FINITE  recognizes  that  the 
substructures  are  linear  and  requests  computation  of  the  condensed  stiff- 
ness matrix  and  equivalent  loads  only  once  during  the  entire  solution. 

A nonlinear  MATERIAL  termed  STEEL  is  defined  for  association  with  finite 
elements  describing  the  vessel.  The  material  model  (yield  function,  flow 
rule,  and  hardening  rule)  selected  from  the  FINITE  material  library  is 
called  VON-MISES.  Details  of  this  model  are  given  in  Appendix  B.  The 
properties  TOLERANCE,  ALPHA,  and  MAXINCR,  specified  by  the  user,  unless 
default  values  are  adequate,  control  numerical  algorithms  within  the 
material  model.  Physical  characteristics  of  the  material  are  specified 
by  a stress-strain  function  selected  from  the  library  (SEGMENTAL  in  this 
problem).  This  particular  function  permits  any  type  of  segmental  stress- 
strain  curve  to  be  input  with  a shortened  form  available  for  elasto- 
plastic  behavior  as  in  this  example.  Elements  in  the  junction  structure 
are  declared  materially  nonlinear  simply  by  specifying  MATERIAL  STEEL 
after  stating  the  element  type. 
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6.5  Deep  Tunnel  in  Rock 

A simple  application  of  the  finite  element  method  to  a nonlinear  j 

geotechnical  problem  involves  excavation  of  a tunnel  in  deep  rock.  | 
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Stresses  around  tunnels  In  deep  rock  were  previously  Investigated  with 
the  finite  element  method  by  Reyes  and  Deere  [60]. 

The  particular  problem  chosen  for  presentation  here  is  a 20  ft. 
diameter  tunnel  excavated  in  rock.  The  Mohr-Coulomb  yield  surface  is 
approximated  by  the  analytically  simpler  Drucker-Prager  yield  criterion 
for  plane  strain  as  discussed  in  Appendix  B.  This  material  model  ade- 
quately represents  the  response  of  rock  for  small  strains. 

The  infinite  medium  surrounding  the  tunnel  is  represented  by  a 
finite  element  mesh  extending  five  tunnel  diameters  away  from  the  tunnel 
face  as  shown  in  Fig.  6.5.1b.  Disturbances  due  to  the  tunnel  are  small 
at  this  distance.  Twenty-eight,  eight-node  quadratic  isoparametric 
elements  were  chosen  to  simulate  the  region  near  the  tunnel.  Ref.  [60] 
employed  480  constant  strain  triangles.  Double  symmetry  planes  permit 
use  of  only  one  quarter  of  the  continuum. 

Loads  applied  to  the  model  simulate  the  construction  sequence.  In 
situ  stresses  of  1000  psi  vertically  and  250  psi  horizontally  are  intro- 
duced by  appropriate  loads  applied  along  the  exterior  boundaries  and  the 
tunnel  face.  Under  these  stresses,  no  elements  are  yielded.  Incremental 
loads  applied  to  the  tunnel  face  gradually  reduce  the  tractions  to  zero 
(representing  completion  of  the  excavation). 

Convergence  was  rapid  with  the  tangent  stiffness  matrix  updated 
before  the  second  iteration  of  each  step.  Most  of  the  FINITE  input  data 
is  given  in  Fig.  6.5.2.  ; 

The  spread  of  plastic  zones  is  indicated  in  Fig.  6.5.1a.  The 
redistribution  of  stresses  due  to  yielding  is  shown  in  Fig.  6.5.3.  For 
both  this  and  the  previous  study,  yielding  was  found  to  spread  at  a 
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forty-five  degree  angle  away  from  the  tunnel  face.  Differences  in  the 
plastic  zones  occurring  In  this  study  and  those  of  Ref.  [60]  are  attributed 
to  the  different  element  meshes.  However,  the  differences  are  negligible  for 
practical  purposes. 

6.6  Thick  Penetrated  Plate 

This  final  example  Illustrates  the  solution  of  a large  (by  current 
standards)  nonlinear  three  dimensional  structure.  The  multiple  restart 
feature  of  the  FINITE  system  was  employed  to  partition  the  total  solution 
process  Into  manageable  segments. 

A 180°  section  of  a thick  circular  aluminum  plate  Is  shown  In  Fig. 

6.6.1.  The  plate  has  a diameter  of  14  In.  and  Is  1.25  In.  thick.  Three 
penetrations  with  chamfered  edges  are  placed  90°  apart.  The  Von  Mises 
material  model  with  an  elastic-plastic  uniaxial  material  stress-strain 
curve  was  selected  to  represent  material  behavior.  A total  of  80, 

20-node  Isoparametric  solid  elements  model  the  plate.  Twenty-eight 
elements  on  the  plate  periphery  were  declared  linear  to  reduce  computation 
time.  Standard  3x3x3  Gaussian  Integration  within  these  elements  Is 
adequate.  A nimierlcal  Integration  order  of  3 x 3 x 5 (5  layers  of  9 
points)  In  the  nonlinear  elements  permits  early  detection  of  yielding 
through  the  depth  caused  by  large  bending  stresses.  The  structure  has 
665  nodes  with  a total  of  1995  degrees  of  freedom.  Nodes  along  the 
outside  edge  of  the  lower  surface  were  completely  restrained  against 
translation.  Input  data  Is  summarized  In  Fig.  6.6.2. 

The  plate  was  analyzed  for  a 100  psi  uniform  pressure  from  which 
the  pressure  at  first  yield  was  computed  as  760  psi.  Three  load  Increments 
were  then  applied  corresponding  to  total  pressures  of  800,  1500,  and 
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2500  pst.  The  elastic  stiffness  was  used  for  the  100,  800,  and  1500 
psi  load  levels.  A new  tangent  stiffness  was  generated  before  the  second 
iteration  of  step  4 (2500  psi). 

The  multiple  restart  feature  of  FINITE  was  utilized  to  solve  each 
load  step  in  a separate  computer  run.  After  the  solution  for  a step 
converged,  the  displacements  and  stresses  were  printed  to  determine  the 
load-deformation  response  and  the  extent  of  yielding.  The  analysis  was 
restarted  via  the  commands  shown  in  Fig.  6.6.3  for  each  subsequent  step. 
The  first  command  requests  POLO  to  execute  the  FINITE  subsystem  with  data 
bases  containing  previously  computed  results  for  the  plate.  The  plate 
structure  is  accessed  and  a new  load  increment  defined.  Convergence 
tests,  tolerances,  output  commands,  etc.  can  also  be  given.  Input  data 
shown  in  the  figure  also  illustrates  the  DESTROY  command  that  allows  the 
analyst  to  release  space  on  secondary  storage  containing  results  no  longer 
required.  Released  space  is  re-used  by  FINITE  to  eliminate  unnecessary 
expansion  of  secondary  storage  requirements.  The  final  commands  request 
computation  of  displacements  and  output  of  results  for  a load  step. 

Segmentation  of  the  total  analysis  into  multiple  computer  runs  is 
mandatory  for  large  structures.  Even  on  large  scientific  computers  the 
exposure  time  (elapsed  wall  clock  time)  can  become  excessive  in  multi- 
processing environments.  The  plate  structure  solution  discussed  in  this 
example  required  over  28  hours  of  elapsed  time  on  a Burrough's  B-6700 
computer.  CPU  times  were  272,  35,  173,  and  374  minutes  for  steps  1 to 
4 respectively. 

The  spread  of  plastic  zones  is  indicated  for  the  various  load  levels 
in  Fig.  6.6.4.  Yielding  begins  between  the  two  cutouts  due  to  large 
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bending  stresses  then  quickly  spreads  over  the  top  and  bottom  surface 
before  penetrating  through  the  depth.  Similar  behavior  Is  observed  In 
plastic  analysis  of  beams.  Extensive  yielding  at  2500  psi  applied  load 
Indicates  the  useful  capacity  of  the  plate  has  been  exceeded. 
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CHAPTER  7 

SUMMARY  AND  CONCLUSIONS 

7.1  Summary 

This  study  examined  the  principles  of  nonlinear  finite  element 
analysis  and  incorporated  them  with  the  principles  of  modern  software 
engineering  to  produce  a readily  usable  tool.  Equal  emphasis  in  the 
presentation  was  placed  on  the  theoretical  formulation  and  the  software 
engineering  aspects  of  finite  element  analysis. 

The  governing  equations  of  nonlinear  continuum  mechanics,  usually 
expressed  in  tensor  notation,  are  presented  in  matrix  form  suitable  for 
direct  application  in  the  finite  element  method.  Both  geometric  and 
material  nonlinearities  are  considered.  The  Lagrangian  formulation  was 
chosen  over  the  Eulerian  formulation  for  reasons  of  computational  effic- 
iency and  convenience  for  the  analyst.  Strain-displacement  relations 
and  stress  definitions  appropriate  for  the  Lagrangian  description  are 
presented  and  physically  interpreted.  Geometric  relationships  derived 
between  the  initial  and  deformed  structural  configuration  enable  the 
analyst  to  transform  between  stress  definitions  as  convenient  during 
analysis. 

Classical  equations  of  the  Lagrangian  formulation  are  cast  into  the 
finite  element  method.  Specific  matrices  given  for  large  deformation  of 
a three  dimensional  continuum  illustrate  the  procedure.  Solution  of  the 
resulting  nonlinear  algebraic  equations  is  accomplished  with  the  Newton- 
Raphson  method  and  its  various  modifications.  A technique  to  account  for 
nonconservative  loads  is  discussed.  Substructuring  and  static  condensation 
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are  proposed  to  reduce  the  cost  of  nonlinear  analysis  by  decreasing  Input 
preparation  effort  and  by  eliminating  consideration  of  linear  regions  of 
the  structure  during  nonlinear  computations. 

The  characteristics  and  features  required  of  a general  nonlinear 
finite  element  computer  system  from  the  analyst's  viewpoint  are  examined. 
These  Include  techniques  for  user-system  communication,  description  of  a 
nonlinear  structural  model  via  substructuring,  analysis  restart,  and  error 
recovery.  Characteristics  of  an  Interface  mechanism  between  the  system 
and  element-material  model  libraries  are  described  that  minimize  the 
effort  to  Implement  new  features. 

An  Implementation  of  these  features  within  the  FINITE  structural 
mechanics  system  Is  described.  The  role  of  the  POLO  supervisor  In 
supporting  problem-orlented-language  translation,  data  management,  and 
dynamic  memory  allocation  are  briefly  examined.  Examples  of  hlerarchal 
data  structures  for  representing  data  associated  with  nonlinear  analysis 
Illustrate  typical  constructs  employed  In  FINITE.  Details  of  the  system- 
library  Interface  scheme  designed  during  this  study  are  presented  to  show 
that  finite  elements  and  material  models  can  be  Isolated  from  the  system 
and  each  other  thereby  producing  easily  modifiable  software.  The  overall 
processes  of  nonlinear  analysis  Including  data  structure  Initialization 
and  solution  for  a single  step  are  described. 

A number  of  nonlinear  example  analyses  demonstrate  features  of  the 
FINITE  system.  User-defined  substructuring,  condensation,  and  sophisti- 
cated data  generators  reduce  analysis  effort  and  cost.  Both  geometric  and 
material  noni Inearl  ties  are  considered  In  the  examples. 
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7.2  Conclusions 

The  basic  equations  of  nonlinear  continuum  mechanics  can  be  expressed 
in  matrix  form  familiar  to  most  structural  engineers  and  can  be  readily 
incorporated  within  the  finite  element  method. 

Proper  design  and  implementation  of  computer  software  for  finite 
element  analysis  has  now  become  recognized  as  an  essential  component  in 
the  overall  analysis  process.  User-oriented,  general  purpose  software 
makes  the  latest  technological  advances  readily  available  to  the  engineer- 
ing profession. 

Engineering  software  support  systems,  commonly  referred  as  super- 
visors, are  a necessary  and  natural  component  of  sophisticated  structural 
analysis  software.  The  use  of  a supervisor  centralizes  data  and  memory 
management  functions  thereby  reducing  development  time  and  cost.  In 
addition,  a supervisor  provides  a convenient  environment  in  which 
developers  may  incorporate  good  software  engineering  principles. 

A comprehensive,  multi-level  substructuring  and  static  condensation 
capability  makes  feasible  the  nonlinear  analysis  of  many  structures 
heretofore  considered  impossible  or  economically  infeasible.  The  full 
potential  of  these  techniques  are  not  realized  unless  multiple  levels  of 
substructuring  with  condensation,  including  loads,  is  supported  auto- 
matically by  the  software. 

The  useful  life  span  of  a nonlinear  finite  element  system  depends 
primarily  on  the  flexibility  designed  into  the  interface  between  the  system 
and  element-material  model  libraries.  A significant  portion  of  the  effort 
required  to  incorporate  new  elements  and  models  can  be  removed  from  the 
developer  and  placed  at  the  system  level. 
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APPENDIX  A 


ELEMENT  TANGENT  STIFFNESSES 


A.l  Formulation 

In  the  Lagrangian  formulation,  forces  exerted  by  a distorted  element 
on  its  nodes  are  given  by 


{IF}  = 


[B]'{o)  dV 


(A.l) 


where  the  [B]  matrix  relates  differential  changes  of  Green  strain  to 
differential  changes  of  element  nodal  displacements,  {a}  is  the  vector 
of  2nd  Piola-Kirchoff  stress  components.  Both  [B]  and  (a)  are  functions 
of  element  nodal  displacements  and  element  local  coordinates  for  general 
nonlinear  analysis.  Integration  is  performed  over  the  initial  configura- 
tion of  the  element.  Components  of  {IF}  are  always  directed  along  the 
undeformed  local  element  axes. 

The  element  tangent  stiffness  matrix,  [Ky],  provides  a linear 
relationship  between  differential  internal  forces  and  nodal  displacements 
such  that 


d{IF}  = [IC^]  d{U} 


(A.2) 


where  {U}  is  the  vector  of  element  nodal  displacement  components.  The 
total  differential  of  {IF}  is  given  by 


d{IF}  = 


d[B]^{a}  dV  + [B]^  d{o}  dV 

V Jv 


(A.3) 


The  second  integral  is  readily  expressed  in  terms  of  displacement  differ- 
entials by  the  following  well  known  substitutions 
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d{or)  = [D^]  d{e}  (A.4) 

die}  = [B]  dlU)  (A. 5) 

which  yield 

[Bf  dlo)  dV  = [f  [B]^  [Dt-]  [B]  dV]  dlUl  (A. 6) 

Jv  Jv  ' 

The  right  integral  is  termed  the  [K®]  portion  of  the  element  tangent 
stiffness.  For  geometric  nonlinear  analysis,  [B]  in  the  above  equation 
is  a function  of  the  nodal  displacements  (U).  For  small  displacement 
analysis,  [B]  is  the  familiar  differential  operator  dependent  on  the 
element  coordinates. 

Computation  of  the  element  stiffness  is  conveniently  split  into 
nodal  stiffness  submatrices  Submatrix  (ij)  relates  differential 

forces  at  node  i to  differential  displacements  at  node  j.  The  matrix 
[k9.]  is  then  given  by 

[K®j]  = f [0^]  [Bj]  dV  (A.7) 

The  three  dimensional  expression  for  [B^]  including  all  nonlinear 
terms  is  given  by  (see  Eq.  3.2.7) 
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(A.8) 


where  x,  y,  and  z subscripts  ii!^)ly  differentiation.  All  derivatives  of  N 
are  for  the  1^^  node  shape  function  only. 

Expansion  of  the  first  integral  In  Eq.  A.3  requires  an  expression 
for  d[B^].  From  the  definition  of  deformed  coordinates  {X'} 


X'  = z (X^  + u^) 

Y*  = z (Y^  + v^)  1 = 1,2....  nn 

Z'  • E (Z^  + w^) 


(A.9) 


their  derivatives  are  easily  obtained.  The  number  of  element  nodes  Is  nn 
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For  example.  X'  is  given  by 


»•  - 3X1 

y ay 


nn  dN 


iiii  Oiij 

■ ^ WT  * "0 


(A.IO) 


Dependence  of  [B^]  on  the  nodal  displacements  becomes  evident  with 
substitution  of  all  derivatives  similar  to  Eq.  A.IO  into  Eq.  A. 8.  The 
[B^]  matrix  is  split  into  two  matrices,  [b!|]  independent  of  nodal  dis- 
placements, and  [B^*"],  dependent  on  nodal  displacements,  by  noting  that 
cross  derivatives  of  nodal  coordinates  are  zero.  For  example. 


r — - X = 0 

1 3y  N 


(A.ll) 


Similarly,  X = Y = Z * 1.  Thus,  in  taking  differentials  of  [B.], 

A V fc  I 


dCB"!]  = 0. 


Differentials  of  each  term  in  [b”'"]  are  obtained  using  the  chain 


rule.  For  example,  is  given  by  • 


<»nl  m du. 

I, I 8u.  j 3u.  'ax  ax  j 

V J 


(A.12) 


Substitution  of 5:(aN|^/ax)  U|^  = au/ax  and  differentiating  with  respect  to 

the  nodal  displacements  yields 

Ml  aN. 

-•SlS  ' h'  ^ sT  “x>  "“j  ' 

J 


aN.  aN. 

^ du. 
ax  ax  J 


(A.13) 


If  5:(aNyax)dUj  is  denoted  du^,  then  dBljl*"^  = (aN^/ax)dUj^. 

NL 

Differentials  of  other  [B^  ] components  are  obtained  in  a similar  way. 
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Differential  Internal  force  components  at  the  1^"  node  can  now  be 
written  more  explicitly.  For  example,  d(IF][)  Is  given  by 


"x  ““x  * "y  "y  % * “z  "z  "“z  * 

"xy  <"x  % * "y  <'“x>  * 

"xz  ("z  ““x  * “x  *'z>  * 

Oyz  ("y  ‘‘“j  ♦ Nj  ] <1* 


(A.M) 


The  differential  force  expressions  must  be  factored  Into  a matrix  product 

to  express  In  stiffness  form.  The  terms  du  , du  , etc.  are  first  written 

X y 

In  matrix  form  as 

d{0r}gxi  = Z [G^]  d{U}j  (A.15) 

where  dllT}^  * [du^,  dv^.  dw^^,  du^,  dv^,  dw^,  du^,  dv^,  dw^]  and 


[G.]  = 3 3y 


(A. 16) 


Is  a 3 X 3 Identity  matrix.  Define  a matrix  [M]  such  that 


[M]  - IgOjjyl  IgCy 


3®xy  1 
1. 

's^xz 

3‘^y  1 
1. 

'3"yz 

3®yz  1 

1 

1 

1 ” 

(A.17) 
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then,  simple  multiplication  shows  that  differential  forces  at  node  1 due 
to  differential  displacements  at  node  j for  the  first  integral  of  Eq.  A.3 
are  given  by 


d{IF^}  * [ 


[M]  [Gj]  dV]  d{Uj} 


(A.18) 


Therefore,  the  total  element  tangent  stiffness  submatrix  (1j)  Is  given 
by 


^[G^f  [M]  [Gj]  + [B^f  [D^]  [B.]  dV 


(A.19) 


A. 2 Truss  Element 

The  three  dimensional  truss  element  with  two  nodes  provides  a simple 
example  to  Illustrate  computation  of  the  element  tangent  stiffness.  Only 
one  local  coordinate,  X,  Is  required.  Displacements  along  the  element  In 
terms  of  nodal  displacements  are 


U(X)  = U2 

V(X)  = + N2  V2  (A.20) 

W(X)  = + N2  W2 


where  the  interpolation  functions  are  defined  by 

(X)  = (L-X)/L 
N2  (X)  = X/L 

The  [6^]  and  [62]  are  then 

[G^]  = -1/L  [I3] 

[Gg]  = 1/L  [I3] 


(A.21) 


(A.22) 
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since  all  derivatives  with  respect  to  Y and  Z vanish.  Similarly,  Is 
the  only  non- zero  stress  component;  thus, 

[M]  » [I3]  (A.23) 

The  matrices  [B^]  and  [Bg]  are  simple  as  Is  the  only  strain  component. 

[B,]  • (^)  [Xi.  y;.  zy  (A.  24) 

Computation  of  the  deformed  coordinates  yields 

X;  = ^ [N^(X^+U^)  + NgCXg+Ug)]  = 1+(U2-U,)/L  (A.25) 

y;  = (V2  - V^)/L 
z;MW2-Wi)/l 

Therefore, 

[B^]  = -1/L  [(1+(U2-Up/L),  (V2-Vp/L,  (W2-Wp/L]  (A.26) 

[B2]  = -[B^] 

For  material  nonlinearity,  [D^]  Is  a 1 x 1 matrix  containing  the 
uniaxial  stress-strain  curve  slope,  Ej.  For  this  element,  the  [G]  and  [B] 
matrices  are  Independent  of  the  element  coordinates,  therefore  the  products 
In  Eq.  A.19  can  be  taken  from  under  the  Integral  sign.  The  resulting 
6x6  tangent  stiffness  matrix  Is  given  below: 


[K^]  = (A  * o^)/L 


l3_  '^-h_ 

"^3  1 ^3 


+ (A  * E^)/L  [B^,B2]' 


(A.27) 
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A. 3 Truss  Element  Implementation  In  FINITE 

To  complete  the  presentation  of  element  tangent  stiffness  generation, 
the  library  input  tables  and  FORTRAN  subroutine  for  the  three  dimensional 
truss  element  are  given.  Fig.  A.l  lists  the  POL  statements  required  to 
define  the  spacetruss  element  to  FINITE.  Most  of  this  data  is  self- 
documenting;  detailed  explanations  are  therefore  omitted. 

The  FORTRAN  subroutine  called  by  FINITE  to  generate  the  tangent 
stiffness  matrix  is  shown  in  Fig.  A. 2.  This  routine  generates  the  tangent 
stiffness  for  a rod,  planetruss,  or  spacetruss  element.  The  labeled 
COMMON  area  in  FINITE  element  routines  minimizes  the  number  of  local 
variables.  Variables  in  this  COMfXlN  space  are  not  passed  between  element 


routines. 
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APPENDIX  B 

PLASTICITY  MATERIAL  MODELS 


B.l  General 

Two  yield  functions  commonly  used  in  plasticity  material  models 
were  implemented  in  FINITE  during  this  work.  The  Von  Mises  yield  criterion 
with  isotropic  hardening  adequately  predicts  the  behavior  of  many  ductile 
metals  subjected  to  monotonic  loading.  The  Drucker-Prager  criterion  [18] 
approximates  the  flohr-Coulomb  yield  function  and  is  relevant  for  some 
cohesive  materials  including  concrete  and  rock.  Both  models  are  standard 
in  most  nonlinear  finite  element  systems  and  are  used  in  this  work  to 
a)  verify  the  material  modeling  processors  of  FINITE,  and  b)  demonstrate 
the  system's  problem  solving  capability.  A brief  review  of  the  formulation 
and  computational  procedures  is  presented  as  well  as  listings  of  the  tables 
and  subprograms  for  incorporating  the  models  into  FINITE.  Complete 
discussions  of  plasticity  theory  as  oriented  toward  finite  element  applica- 
tions can  be  found  in  Refs.  [11,  48,  49,  50,  74]. 

The  implementation  of  these  models  in  FINITE  is  readily  adaptable 
to  other  yield  criterion  and  associative  flow  rules  that  follow  general 
plasticity  theory.  Special  numerical  techniques  used  here  include  the 
"subincrement"  method  that  prevents  artificial  hardening  due  to  replacement 
of  differential  quantities  by  finite  increments.  Additional  corrective 
techniques  are  employed  during  the  transition  between  elastic  and  plastic 
states. 

B.2  Formulation 

Experimental  results  have  verified  the  negligible  effect  of  hydro- 


static stress  on  yielding  in  ductile  metals.  Von  Mises  proposed  that 


yielding  of  a material  under  a complex  state  of  stress  is  linked  to  the 
distortional  strain  energy  density,  i.e.,  energy  due  to  change  of  shape  not 
change  of  volume.  At  initiation  of  yielding,  the  distortional  strain 
energy  under  complex  stress  states  is  assumed  equal  to  that  for  simple 
shear.  Equating  the  two  expressions  for  strain  energy  defines  a yield 
function  given  by 

^ - o (B.l) 


where. 


and 


J = 1 (S^  + + S^)  + 

2 2 X y 2 xy  xz  yz 


(B.2) 


(B.3) 


= (a^  °y 

a is  the  yield  stress  of  the  material  in  uniaxial  tension.  Expressed  in 
terms  of  principal  stresses,  defines  a cylindrical  surface  in  stress 
space  with  the  generator  along  the  line  0-1=02  = and  radius  of  length 
o /2/3.  The  intersection  of  the  cylinder  with  the  0^02  plane  defines  the 
well  known  Von  Mises  ellipse  for  two  dimensional  states  of  stress.  These 
surfaces  are  shown  in  Fig.  B.l. 

The  Drucker-Prager  yield  function  is  similar  to  the  Von  Mises  criterion 
but  incorporates  the  effect  of  hydrostatic  stress.  The  yield  function  is 
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u 

I ; 

I 

!.  , 

expressed  by 

fj)p  " '^rn  ‘^2  " 

where  a and  K are  material  constants  dependent  upon  the  cohesion  strength, 
C,  obtained  from  an  unconfined  compression  test  and  « , the  angle  of 
Internal  friction.  The  Drucker-Prager  yield  function  plots  as  a right 
circular  cone  with  the  same  generator  axis  as  the  Von  Mises  cylinder  and 
approximates  the  Mohr-Coulomb  pyramid  yield  surface.  The  cone  diameter  Is 
proportional  to  the  hydrostatic  compression. 

Three  sets  of  coefficients  a and  K have  been  proposed  [60].  The 
first  set  Is  obtained  by  forcing  the  cone  to  pass  through  edges  of  the 
Mohr-Coulomb  pyramid  determined  from  a triaxial  test.  The  second 
coefficients  approximate  the  pyramid  for  the  axial  extension  test.  The 
coefficients  for  these  two  sets  are 

_ 2 SIN  ♦ 

O = 

' ^ (3-SIN  $) 

= 2 SIN  $ 

^ /J  (3+SIN  ♦) 

Drucker  and  Prager  gave  the  following  constants  for  the  special  case  of 
plane  strain 

TAN  ♦ „ _ 3 C /„ 

Oo 5 T75-  Ko  “ = — 5 YJJ  vB.o; 

(9  + 12TAfr  ^ (9  + 12TAN^ 

For  given  values  of  C and  the  constants  of  B.8  have  the  least  values  of 
a and  K and  therefore  are  the  more  conservative.  The  analyst  may  choose 
any  set  of  coefficients  through  the  material  model  property  COEFFICIENTS. 
Values  1,  2,  and  3 for  this  property  correspond  to  coefficients  defined 
by  Eqs.  B.6,  B.7,  and  B.8  respectively. 

i 

I 


6 C COS  » 
^ (3-SIN  «) 


(B.6) 


= 


6 C COS  $ 


^ vT  (3+SIN  ♦) 


(B.7) 


All  possible  elastic  stress  states  lie  inside  the  cylinder  and  cone 
(F  < 0).  Stress  states  on  the  yield  surface  (F  = 0)  correspond  to  yield- 
inQ  or  a plastic  condition.  Materials  cannot  inaintain  stress  states 

outside  the  yield  surface  (F  > 0). 

Once  the  stress  state  reaches  the  yield  surface  during  incremental 
loading,  a flow  rule  is  required  to  relate  subsequent  changes  in  stress 
to  increments  of  strain.  Most  theories  assume  a linear  relationship 
between  infinitesimal  stress  and  strain  increments.  By  examining  the 
implications  of  work  hardening  and  the  ideally  plastic  material,  Drucker 
concluded  that  plastic  strain  increments  are  proportional  to  the  gradient 
of  the  yield  surface.  Thus, 

{dcp}  = dx  (al  (B.9) 

where  {a}  contains  components  of  the  outward  normal  to  the  yield  surface 
in  stress  space  and  dx  is  a positive  constant.  This  relation  is  termed 
the  normality  principle. 

A hardening  rule  is  required  to  predict  changes  in  the  yield  surface 
due  to  plastic  deformation  of  the  material.  As  work  hardening  enables 
the  material  to  maintain  higher  states  of  stress,  the  yield  surface  must 
be  modified  accordingly  to  keep  the  stress  state  on  the  surface.  The 
simplest  form  of  hardening  assumes  that  the  yield  surface  retains  its 
original  shape  and  expands  uniformly  to  reflect  the  new  yield  stress. 

This  type  hardening  rule  is  termed  isotropic  and  ignores  the  Bauschinger 
effect.  Other  hardening  rules  distort  the  yield  cylinder  into  a cone,  for 
example,  or  rotate  and  translate  the  cylinder  as  a rigid  body  in  stress 
space.  For  the  Von  Mises  model,  isotropic  hardening  is  achieved  by 


i 


PU 


equating  the  plastic  work  done  by  actual  stresses  to  the  plastic  work  done 
by  a specimen  of  the  material  under  uniaxial  tension.  No  strain  hardening 
Is  permitted  for  the  Drucker-Prager  model  Implemented  here. 

With  the  yield  function,  flow  and  hardening  rules  established, 
derivation  of  the  constitutive  relation  valid  In  the  plastic  range  Is 
relatively  straightforward.  For  a given  Infinitesimal  Increment  of  strain, 
{de},  composed  of  elastic,  (de  },  and  plastic,  (de  },  portions  the  follow- 

6 P 

Ing  relations  hold 


(do)  = [D]  {de^} 


{do}  = [D]  {del  - [D]  {dCp} 


(B.IO) 


(B.ll) 


where  [D]  Is  usual  matrix  of  elastic  constants  for  Isotropic  materials 
and  {dcp}  Is  computed  by  Eq.  B.9.  When  stresses  are  such  that  F = 0, 
the  differential  must  also  be  zero.  In  the  Von  Hises  criterion,  for 


example. 


F = dF  = {da}  - dS’  = 0 
3a 


da  = {|^}  (da) 


(B.12) 


(B.13) 


where  {3F/3a}  Is  the  stress  gradient  normal  to  the  yield  surface.  Thus, 
(a)  of  Eq.  B.9  Is  given  by  {3F/3a}  . Equating  the  plastic  work  done  by 
the  stress  state  causing  yield  to  the  plastic  work  done  by  the  uniaxial 
yield  stress  shows  that 


{a}  {dcp}  » a dcp 


(B.14) 


Substitution  of  the  plastic  strain  expression  from  the  flow  rule,  Eq.  B.9, 
yields 

dx  {o}^  {a}  = a dcp  (B.15) 

Euler's  theorem  for  homogeneous  functions  shows  that  {a}^{a}  = 0 from 
which 


dx  = dcp  (B.16) 

Pre-multi plying  Eq.  B.ll  by  {a}^  gives 

{a}^  {do}  = do  = 

(a)^  [D]  {dc}  -{a}^  [D]  dl  {a}  (B.17) 

The  increment  of  uniaxial  stress,  do,  is  given  by 

do  = H’  dip  (B.18) 

is  the  slope  of  the  uniaxial  tension-plastic  strain  curve 
from  test  results.  Substituting  into  Eq.  B.17  and  solving  for 

dx  = dl  = {de}  (B.19) 

P H'  + {a}'  [D]  {a} 

Substitution  of  this  expression  into  Eq.  B.ll  provides  the  desired  relation 
between  stress  and  strain  in  the  plastic  range 


where  H' 
obtained 


133 


u 

u 

LJ 


i 

L 


n 


(_j 


The  above  expressions  are  used  to  evaluate  the  change  In  stress 
state  for  a finite  Increment  of  strain  {Ae}  as  follows 


{La) 


R [D]  {Lt) + 


tLz 


[Or]  {de> 
RAe  ' 


(B.21) 


where  R Is  the  portion  of  the  strain  Increment  required  to  reach  the 
yield  surface  from  a previously  elastic  stress  state.  If  the  previous 
stress  state  was  on  the  yield  surface,  R = 0.  For  small  Increments  of 
load  resulting  in  small  strain  Increments,  the  above  expression  Is 
approximately  equal  to 


{Aa}  = R [D]  (Ae)  + (1-R)  [D^]  {Ae}  (B.22) 


This  approximation  Is  Inadequate  for  large  strain  Increments  and  causes 
drifting  from  the  yield  surface  thereby  Introducing  artificial  hardening 
as  shown  In  Fig.  B.3.  Special  numerical  techniques  to  reduce  drifting 
effects  are  described  In  the  next  section. 


B.3  Numerical  Refinements 

The  R Factor  for  transitioning  between  elastic  and  plastic  stress 
states  Is  computed  more  accurately  by  the  process  outlined  below  applicable 
for  all  convex  yield  surfaces.  This  procedure  was  originally  developed 
by  Nayak  [50].  Let  {oq}  denote  an  Initial  stress  state  Inside  the  yield 
surface  as  shown  In  Fig.  B.3.  If  the  entire  strain  Increment  Is  considered 
elastic,  the  new  stress  state  Is  outside  the  yield  surface.  This  condition 
Is  expressed  by 

Fdap})  = Fq  < 0 

Uog}  = [D]  {Ae}  (B.23) 

F({Oq1  +{AOg})  = F^  > 0 
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The  R factor  is  chosen  such  that 


F({oq  + R{aOg})  = 0 


(B.24) 


Simple  linear  interpolation  between  and  F^  provides  the  first 


estimate  for  R 


R,  = -F,/(F,  - F„) 


(B.25) 

But  because  F is  a nonlinear  function  of  the  stress,  the  correction  is 
not  exact. 


F({aQ}  + R^{AOg})  = Fg  M (B.26) 

A small  chanQO  in  R-j  is  necessary  to  produce  a chanQe  in  F,  i.e.,  F uF  = 
0.  For  an  instantaneous  position  of  the  yield  surface,  dF  is  given  by 
{a}^{do}  with  a evaluated  at  the  stress  state  {oq}  + Ri{Aag}.  Approxi- 
mating {da}  by  AR^CAOg},  the  change  in  R^  is  given  by 


AR^  = -Fg/CCal'lAOg}) 


(B.27) 


with  the  improved  value  for  R given  by  R = R-j  + aR^.  This  process  has  a 
simple  geometric  interpretation  as  shown  in  Fig.  B.3.  The  projection  of 
the  estimated  stress  state  to  cause  yield,  computed  using  R^ , onto  the 
yield  surface  normal  is  forced  to  zero  by  the  correction  dF. 

Drift  from  the  yield  surface  caused  by  approximating  the  integral 
of  Eq.  B.21  by  the  expression  in  Eq.  B.22  is  minimized  with  the  "subincrement" 
procedure.  The  technique  is  based  on  the  geometrical  interpretation  of 


the  [Dj]  matrix  as  shown  in  Fig.  B.3.  This  matrix  forces  the  new  stress 


state  to  remain  on  the  yield  surface  for  infinitesimal  strain  increments. 


I 


^ --^1*1***  ■'  - v . 


When  finite,  rather  than  Infinitesimal,  strain  Increments,  occur  during 
actual  problem  solution  considerable  drift  from  the  yield  surface  may 
result.  To  minimize  drift,  the  portion  of  the  strain  Increment  not  taken 
elastically,  (1-R)  (Ac),  Is  subdivided  Into  M equal  Increments.  The 
Integral  In  Eq.  B.21  Is  then  approximated  by 


rAe  M 

[D^]  {dc}  = (1-R)  I [ap]{Ae}/K 

RAc  ^ 


(B.28) 


In  which  [Dy]  Is  updated  at  the  end  of  each  subincrement.  The  number  of 
subincrements  Is  based  upon  the  magnitude  of  the  Initial  deviation  from 
the  yield  surface  for  the  entire  strain  Increment  taken  elastically  or 
the  deviation  computed  using  the  [D^]  Incorporated  In  the  current  element 
stiffness  matrix  If  the  material  has  previously  yielded.  If  the  Initial 
deviation  Is  denoted  F^,  then  the  number  of  Increments  Is  given  by 


M = F,/(oo)  < (Von  Mises) 

1 lUoX 


M = F^/(aK)  < (Drucker-Prager) 


(B.29) 


where  a Is  a fraction  of  the  current  yield  stress.  For  small  deviations 
only  a few  subincrements  are  required.  A limit  on  the  maximum  number  of 
Increments  Is  necessary  to  detect  Instability  of  the  procedure  (a  large 
number  of  subincrements  Indicates  overall  analysis  divergence). 

Three  user  supplied  material  model  properties  control  the  subincre- 
ment procedure.  The  property  ALPHA  corresponds  to  o In  Eq.  B.29  and  has 
a default  value  of  0.05.  The  maximum  number  of  Increments  Is  given  by  the 
property  MAXINCR  and  has  a default  value  of  30.  DIVERGE  Is  the  number  of 


increments  above  which  the  solution  is  said  to  be  diverging  (default  value 
of  500).  If  the  number  of  subincrements  is  greater  than  HAXINCR  but  less 
than  DIVERGE,  the  number  of  increments  used  is  MAXINCR.  One  final 
property,  TOLERANCE,  is  used  to  determine  if  a state  of  stress  is  on  the 
yield  surface  F.  If  the  yield  function  F has  an  absolute  value  less  than 
TOLERATiCE  * o for  the  Von  Mises  criterion  or  TOLERANCE  * K for  the  Drucker- 
Prager  criterion,  then  the  state  of  stress  is  assumed  to  lie  on  the  yield 
surface.  The  default  value  of  TOLERANCE  is  0.001. 

Even  with  the  subincrement  process,  updated  stresses  at  the  end  of 
each  subincrement  will  not  satisfy  the  yield  criterion.  Let  {o'}  denote 


updated  stresses  following  each  subincrement,  then 


(B.30) 


A final  corrective  stress  (Ao')  computed  by  a procedure  similar  to  that  in 


(B.31) 


where  (a)  is  evaluated  for  the  stresses  {o'}. 


B.4  Computational  Procedure 

The  major  steps  for  updating  stresses  given  an  increment  of  strain 
are  summarized  below. 

1)  Subtract  that  portion  of  the  incremental  strain  vector  due  to 
initial  strain  effects  to  obtain  the  increment  of  mechanical 
strain,  {Ae}.  Using  the  constitutive  matrix  for  the  point 
incorporated  in  the  current  element  stiffness,  either  [0]  or 
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[Dy],  compute  the  trial  stress  increment  (AOg}.  Add  these  to 
the  previous  total  stresses  to  obtain  a new  trial  stress 
{o'}  = {oq}  + {AOg}. 

2)  Evaluate  the  yield  function  F({o'})  * F^.  If  the  previous  state 
of  stress  was  plastic,  F({oq})  = Fq  * 0,  set  R = 0 and  go  to 
step  4.  If  Fq  < 0,  the  new  stress  state  is  still  elastic. 

Trial  stresses  {o'}  are  actual  stresses.  Omit  remaining  steps. 

If  F^  >0  and  Fq  < 0,  the  R factor  is  required.  Go  to  the  next 
step. 

3)  Compute  R^  = to  {oq}  and  re-evaluate 

the  yield  function  to  obtain  Fg.  Determine  {a}  and  compute 

R = R^-Fg/Ial^lAOg}.  Generate  the  new  trial  stress  as 
{o'}  * {o-}  + R{Ao  }.  The  yield  function  should  be  very  near 

U C 

zero  for  these  stresses. 

4)  Compute  M,  the  number  of  subincrements.  M = F^/(oo)  + 1 or 
M * F^/(aK)  + 1.  Compute  the  subincrement  strain  increment, 

{Ae"}  ={Ae}/M,  and  the  corresponding  elastic  stress  increment, 
{Ao"}  = [D]  {Ae"}.  Repeat  steps  5-7,  M times. 

5)  Invoke  the  material  stress-strain  function  to  return  H'  for  the 
total  accumulated  uniaxial  plastic  strain.  Omit  this  step  for 
the  Drucker-Prager  model. 

6)  For  the  current  stresses,  {o'},  compute  {a}  and  the  plastic 
multiplier  d using  Eq.  B.19.  If  dx  < 0,  the  material  is  unloading 
elastically  from  a plastic  state.  Compute  the  new  stress 

using  the  [0]  matrix  and  omit  remaining  subincrements. 
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I For  positive  dx,  compute  the  plastic  strain  increments, 

{ACp},  and  adjust  the  elastic  increment  {Ao")  to  reflect  them. 

; {o'}  = {o'}  + {Ao"}  - dx[D]  {a}. 

7)  Update  the  uniaxial  plastic  strain  to  include  dx.  If  H'  = 0, 
adjust  stresses  using  Eq.  B.31.  For  a strain  hardening  material, 

! compute  the  new  yield  stress  o by  evaluating  F({o'})  with 

j previous  o = 0. 

j 8)  The  final  new  stress  state  for  the  material  is  that  computed  for 

the  last  subincrement.  Update  the  material  history  parameters 

i ■ 

■ to  indicate  the  state  is  plastic  and  save  the  accumulated  uniaxial 

! i 

^ plastic  strain. 

B.5  Implementation  in  FINITE 

The  POL  statements  to  define  the  Von  Mises  material  model  and  the 
I SEGMENTAL  stress-strain  function  are  shown  in  Figs.  B.4  and  B.5.  The  two 

corresponding  FORTRAN  subroutines  are  listed  in  Figs.  B.6  and  B.7. 
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Fig.  5.2.1.  Typical  Nonlinear  Software  Organization 


Processing 
Module  n 


Fig.  5.2.2.  Nonlinear  Finite  Element  Software  Organized 
Around  a Data  Base  Manager 


STRUCTURES 


SNAME 

EPSVECPOS 


m. 


Fig.  5.3.1.  Typical  Hlerarchlal  Data  Structure 


POLO  FINITE  Subsystsms 


Fig.  5.4.2.  Data  Structure  for  Interfacing 
FINITE  Subsystems 
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Fig. 


5.6.2.  Initialization  of  Nonlinear  Data  Structure 
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Crown—Node  97 


97  Nodes 
264  Members 


L - Latitudinal 
M- Meridional 
0-  Diagonal 


Span  150  ft 
Rise  15  ft 

E * 29,000  k si 
A = 5 ia* 

All  Members 


(b)  Elevation 


Fig-  6,2.1.  Schwedler  Dome 
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(a)  Substructure  BOTTOM.RING 
(Linear) 
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(c)  structure  DOME 
(Nonlinear) 


I 

Fig.  6.2.2.  Schwedler  Dome  - Substructured  Model 
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Fig.  6.2.3.  FINITE  Input  Data  for  Standard 
Model  of  Schwedler 
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LINE  DC 

DISPLACEMENT, 

STANDARD  MODEL 

IN  INCHES 

SUBSTRUCTURED  MODEL 

Number 

(1) 

U 

(2) 

V 

(3) 

W 

(4) 

U 

(5) 

V 

(6) 

W 

(7) 

2 

.120 

-.392 

-1.301 

-.117 

-.378 

-1.248 

3 

.218 

-.395 

-1.501 

-.213 

-.395 

-1.500 

4 

.246 

-.516 

-2.506 

-.243 

-.508 

-2.452 

5 (Crown)- 

.134 

-.349 

0.464 

-.131 

-.345 

0.458 

LINE  DB 

2 

.038 

-.094 

-.107 

-.050 

-.092 

-.146 

3 

.164 

-.172 

-.587 

-.154 

-.167 

-.542 

4 

.208 

-.290 

-.823 

-.204 

-.288 

-.809 

Fig.  6.2.6.  Schwedler  Dome  — Partial  Loading  Results  for  29  PSF 


Fig.  6.3.3.  Lateral  Deflections  for  P-Delta  Analysis 
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P- DELTA  ANALYSIS  OF  SIMPLE  PLANE  FRAME  STRUCTURE 


t I ELEMENT  BEAM  $ IGNORE  AXIAL  DEFORMATION  IN  BEAMS.  I.E..  AX  - 500. 

TYPE  PLANEFPJUC  E 29600.  TABLE  U N21X62  AX  500. 

COORDINATES 

i 2 432 

! LOADING  FLOOR  ‘ 1.0  KIPS  PER  INCH  ' 

\ ELEMENT  LOADS 

!•  UNIFORM  FORCE  Y W -1.0 

I-  C 

ELEMENT  COLUWl 

TYPE  PLANEFRAME  LARGE  DISPLACEICNTS  E 29600.  TABLE  W W10X49 
COORDINATES 
2 120 
C 

STRUCTURE  FRAME 
NUMBER  OF  NODES  16  ELEMENTS  21 
ELEf^NTS 

1-12  TYPE  COLUMN  ROTATION  Z 90. 

13-21  TYPE  BEAM  ROTATION  SUPPRESSED 
INCIDENCES 

GENERATE  3 IN  X 4 IN  Y FROM  1 5 ADD  4 IN  X 1 IN  Y 
GENERATE  3 IN  X 3 IN  Y AS  13-21  FROM  5 6 ADD  1 IN  X 4 IN  Y 
CONSTRAINTS 
1-4  U V • 0.0 

LOADING  WORKING  'BEAMS  2.8  K/F  PLUS  WIND' 

EXTERNAL  ELEICNT  LOADS 
13-21  FLOOR  0.23333 
NODAL  LOADS 
5 9 FORCE  X 5.0 
13  FORCE  X 2.5 

LOADING  COLUm-LOADS  'EQUIVALENT  COLUMN  LOADS  FROM  GIRDERS' 
NODAL  LOADS 

5 9 13  8 12  16  FORCE  Y -50.4 

6 7 10  11  14  15  FORCE  Y -100.8 
LOADING  P-DELTA  'LOADING  FOR  NONLINEAR  ANALYSIS' 

NONLINEAR 

STEP  1 'WORKING  LOADS'  WORKING  1.0 
STEP  2 'WORKING  + 0.7*AXIAL  'COLUMN-LOADS  0.7 
C 

CONVERGENCE  TEST  NORM  RESIDUAL  LOADS  TOLERANCE  0.1 
TRACE  NONLINEAR  SOLUTION 
MAXIMUM  ITERATIONS  10 
CONTINUE  IF  NONCONVERGENT 
PRINT  RESIDUAL  LOADS 
UPDATE  STIFFNESS  EVERY  2 ITERATIONS 
C 

COMPUTE  NONLINEAR  DISPLACEMENTS  FOR  STRUCTURE  FRAfC, 

LOADING  P-DELTA  STEPS  1 2 

OUTPUT  WIDE  NONLINEAR  DISPLACEMENTS  STRESSES  FOR  STRUCTURE  FRAME, 
LOADING  P-DELTA  STEPS  1 2 

STOP 


Fig.  6.3.4. 


FINITE  Input  Data  for  P-Delta  Analysis 


250  3 0.0201  0.615  258 

2 0.0155  0.548  257 


Average  = 258 


E 29,120  kti 
V 0.3 

o-y  40.54  ksi 


Nonlinear  Region 


Condensed  Cylinder 
Substructure 


139  Nodes 


36  Elements 


8 Node  Quadratic 


Condensed 
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c 

C JUNCTION  OF  SPHERE  AND  CYLINDER. 

C 

STRUCTURE  JUNCTION 

NUMBER  OF  NODES  139  ELEMENTS  38  COORDINATE  POINTS  141 
ELEMENTS 

1-36  TYPE  AXIQZDISOP  MATERIAL  STEEL  MIX  2 HIY  2 

37  TYPE  CYLINDER  ROTATION  SUPPRESSED 

38  TYPE  SPHERE  ROTATION  SUPPRESSED 
INCIDENCES 

168314572  140  141 

2 11  13  8 6 9 10  12  7 140  141 

3 17  19  13  11  15  16  18  12  140  141 


COORDINATES 

ORIGIN  X 2.8125  Y 8.2329 
1 0.  1.5 


GENERATE  2 IN  X 2 IN  Y USING  PATTERN  17-19  Y 22  0 23, 
Y 25-27 


LOADING  UNIT-PRESS  'INTERNAL  PRESSURE  100  PSI' 

ELErtNT  LOADS  FOR  TYPE  AXiq2DIS0P 
1 2 3 5 7 9 11  13  15  UNIFORM  NORMAL  W 0.1  FACE  2 
24-36  BY  2 UNIFORM  NORt^L  W 0.1  FACE  4 
EXTERNAL  ELEMENT  LOADS 
37  38  UNIT-PRESSURE  0.1 

LOADING  PRESSURE  'NONLINEAR  LOADING  CONDITION' 

NONLINEAR 

STEP  1 '700  PSI  INTERNAL  PRESSURE'  COf«INE  UNIT-PRESS  7. 

STEP  2 '800  PSI  INTERNAL  PRESSURE'  UNIT-PRESS  1.0 

STEP  3 '900  PSI  INTERNAL  PRESSURE'  UNIT-PRESS  1.0 

STEP  4 '1000  PSI  INTERNAL  PRESSURE'  UNIT-PRESS  1.0 
STEP  5 '1100  PSI  INTERNAL  PRESSURE'  UNIT-PRESSURE  1.0 
C 
C 

MAXIMUM  ITERATIONS  1 

CONVERGENCE  TEST  NORM  RESIDUAL  LOADS  TOLERANCE  1.0 
TERMINATE  IF  NONCONVERGENT 
UPDATE  TANGENT  STIFFNESS  AT  ITERATIftN  2 
TRACE  NONLINEAR  SOLUTION 
C 

C f 

COMPUTE  NONLINEAR  DISPLACEMENTS  FOR  STRUCTURE  JUNCTION, 

FOR  LOADING  PRESSURE  STEP  1 

OUTPUT  WIDE  NONLINEAR  DISPLACEMENTS  ALL  STRESSES  1-36, 

FOR  STRUCTURE  JUNCTION  LOADING  PRESSURE  STEP  1 f 

STOP 


Fig.  6.4.3.  continued 


Fig.  6.5.2.  FINITE  Input  Data  for  Excavation  of  Deep 
Tunnel  in  Rock 


Redistribution  of  Stresses  in  Deep  Tunnel 
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Fig.  6.6.2.  FINITE  Input  Data  for 

Penetrated  Plate  Structure 


oooooooo 


178 


*RUII  FINITE  FILE  = 20,  , 22 

RESTART  ANALYSIS  OF  THICK  PENETRATED  PLATE. 

DEFINE  THE  3 RD  LOAD  STEP,  INQUEST  TRACE  OF  THE 
NONLINEAR  SOLUTION  PROCESS,  RESET  THE  MAXIHUH  NUMBER 
OF  ITERATIONS  AND  CONVERGENCE  TEST,  THEN  REQUEST 
SOLUTION  FOR  DISPLACEMENTS. 

STRESSES  ARE  REQUESTED  ONLY  FOR  THE  NONLINEAR  ELEMENTS. 

ACCESS  STRUCTURE  PLATE 
C 

LOAOING  PRESSURE 

STEP  3 'TOTAL  PRESSURE  1500  PSI'  UNIT-PRESSURE  700. 

C 

TRACE  NONLINEAR  SOLUTION 
I’AXIMUM  ITERATIONS  10 

CONVERGENCE  TEST  NORM  RESIDUAL  LOADS  TOLERANCE  1.0 
C 

DESTROY  NONLINEAR  RESULTS  FOR  STRUCTURE  PLATE  LOAD  PRESSURE  STEP  1 
COMPUTE  NONLINEAR  DISPLACEMENTS  FOR  STRUCTURE  PLATE, 

FOR  LOADING  PRESSURE  STEP  3 

OUTPUT  WIDE  NONLINEAR  DISPLACEMENTS  FOR  STRUCTURE  PLATE, 

LOADING  PRESSURE  STEP  3 
OUTPUT  WIDE  NONLINEAR  STRESSES  1-2, 

7-9  14-17  21-27  31-38  41-47  50-57, 

60-64  67-71  74-76  FOR  STRUCTURE  PLATE, 

LOADING  PRESSURE  STEP  3 

STOP 


Fig.  6.6.3.  FINITE  Input  Data  for  Penetrated  Plate 
Analysis  Restart 


Uniform  Pressure,  psi 


(a)  Yielding  At  Top  Layer  Integration  Points  For 
Values  of  Applied  Uniform  Pressure 


(b)  Lateral  Deflection  At  Point  A 

Fig.  6.6.4.  Results  for  Thick  Penetrated  Plate 


o n o n o o 


POL  TABLE  DATA  TO  DEFINE  ELEMENT  SPACETRUSS 


DEFINE  ELEMENT  SPACETRUSS 
NUMBER  OF  NODES  2 
NODAL  FREEDOMS 
ALL  U V W 
SUBROUTINE  1 
PROPERTIES 

E REAL  1 . 0 

G REAL  0.0 

NU  REAL  0 . 3 

ALPHA  REAL  0.0 

MASS  REAL  0.0 
DENSITY  REAL  0.0 
DEBUG  LOGICAL  FALSE 
AX  REAL  0 . 0 
END  OF  PROPERTIES 
NUMBER  OF  STRAIN  POINTS  2 
NUMBER  OF  REAL  STRAINS  1 
WORKSPACE  7 WORDS 
NONLINEARITY  GEOMETRIC  MATERIAL 
STRAIN  DECLARATIONS 

ALL  EPX  - SIGMA_X,  DELTAL  - FORCE_X 
LOADS 

TYPE  CONCENTRATED 

P INTENSITY  REAL  -1.0 
L REAL  0.0 

FRACTIONAL  LOGICAL  FALSE 
TYPE  UNIFORM 


END  OF  LOADS 
END  OF  ELEMENT 


n n n n n 


30  CONTINUE 

L2  = L * L 

B(l)  = - { DU/L  + 1.  ) / L 
B(2)  = -DV  / L2 
B(3)  = -DW  / L2 
B(4)  = -B(l) 

B(5)  = -B(2) 

B(6)  = -B(3) 

FACTOR  = AX  * E * L 
DO  50  I = 1,  6 
DO  40  J = 1,  6 

STIFF(I,J)  = STIFP{I,J)  + ( B(I)  * B(J)  * FACTOR  ) 
40  CONTINUE 
50  CONTINUE 

IF  ( .NOT.  DEBUG  ) GO  TO  100 
CALL  ASEERR(  lOUT  ) 
lERR  = 0 

WRITE (lOUT, 1000)  MATNON,  GEONON,  SIGMAX,  DU,  DV,  DW, 
1000  FORMAT (//,3X,13HTRUSS  ELEMENT  ,/ , 3X , 2L5 , 5F10 . 5 ) 
WRITE(IOUT,1020)  E,  AX 
1020  FORMAT (3X,2F20. 6) 

WRITE (lOUT, 1010)  ( (STIFF ( I, J) ,J=1,6) ,1=1,6) 

1010  FORMAT (/,3X,16HSTIFFNESS  MATRIX , 6 (/ , 5X , 6F10 . 3 ) ) 


PASS  REQUESTED  K SUBMATRIX  BACK  TO  FINITE 


100  CONTINUE 

KK  = ( IROW  - 1 ) * 3 
LL  = { JCOL  - 1 ) * 3 

DO  120  I = 1,  M 
DO  110  J = 1,  M 
K(I,J)  = STIFF (KK+I ,J+LL) 
110  CONTINUE 
120  CONTINUE 
C 

RETURN 

END 


Fig.  A. 2. 


continued 


Updated  Yield 
Surface 


Initial  Yield  Surface 


Fig.  B.2.  Normality  Principle  and  Isotropic 
Hardening 
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Fig.  B.3.  Drifting  Errors  and  Transition  Between 
Elastic  and  Plastic  Regions 


POL  TABLE  DATA  TO  DEFINE 


STRESS-STRAIN  FUNCTION  SEGMENTAL 


C ELASTO  PLASTIC  OR  STRAIGHT  LINE  SEGEMENT  APPROXIMATION 


DEFINE  MATERIAL  STRESS_STRAIN  FUNCTION  SEGMENTAL 
SUBROUTINE  1 
PROPERTIES 

YMODULUS  REAL  0 . 0 
NU  REAL  0.0 

SEGMENTS  LOGICAL  FALSE 
NUMPOINTS  INTEGER  0 
STRAINS  VECTOR  REAL  0. 

STRESSES  VECTOR  REAL  0. 

ELASTO_PLASTIC  LOGICAL  FALSE 
COMPYIELD  REAL  0.0 
TENSYIELD  REAL  0.0 
DEBUG  LOGICAL  FALSE 
END  OF  PROPERTIES 
END  OF  FUNCTION 


Example  Stress-Strain  Function 
Definition  Tables  for  FINITE 


Fig.  B.6.  Material  Model  Von-Mises  Subprogram 
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Fig.  B.6.  continued 
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PART  1 - Government 

Administrative  and  liaison  Activities  Navy  (Con't.) 

Naval  Research  Laboratory 
Washington,  DC  20375 
Attn : Code  8400 
8410 
8430 
•8440 
. 6300 
6390 
6380 


Office  of  Naval  Research 
Department  of  the  Navy 
Arlington,  VA  22217 
Attn:  Code  474  (2) 
Code  471 
Code  200 


Director 

Office  of  Naval  Research 
Branch  Office 
666  Summer  Street 
Boston,  MA  02210 


David  W.  Taylor  Naval  Ship  Research 
and  Development  Center 
Annapolis,  MD  21402 
Attn:  Code  2740 
28 
281 


D1 rector 

Office  of  Naval  Research 
Branch  Office 
536  South  Clark  Street 
Chicago,  IL  60605 


U.S.  Naval  Weapons  Center 
China  Lake,  CA  93555 
Attn:  Code  4062 
4520 


Director 

Office  of  Naval  Research 
New  York  Area  Office 
715  Broadway  - 5th  Floor 
New  York,  NY  10003 


Commanding  Officer 

U.S.  Naval  Civil  Engineering  Laboratory 
Code  L31 

Port  Hueneme,  CA  93041 


Director 

Office  of  Naval  Research 
Branch  Office 
1030  East  Green  Street- 
Pasadena,  CA  91106 

Naval  Research  Laboratory  (6) 

Code  2627 

Washington,  DC  20375 

Defense  Documentation  Center  (12) 
Cameron  Station 
Alexandria,  VA  22314 


Naval  Surface  Weapons  Center 
White  Oak 

Silver  Spring,  MD  20910 
Attn:  Code  R>10 
G-402 
K-82 


Technical  Director 
Naval  Ocean  Systems  Center 
San  Diego,  CA  92152 

Supervisor  of  Shipbuilding 
U.S.  Navy 

Newport  News,  VA  23607 

U.S.  Navy  Underwater  Sound 
Reference  Division 
Naval  Research  Laboratory 
P.O.  Box  8337 
Orlando,  FL  32806 


Nayy. 

Undersea  Explosion  Research  Division 
Naval  Ship  Research  and  Development 
Center 

Norfolk  Naval  Shipyard 
Portsmouth,  VA  23709 

Attn:  Dr.  E.  Palmer,  Code  177 


Navy  (Con't.) 


Chief  of  Naval  0|>rrations  C 

Department  of  the  Navy  D 

Washington,  DC  20350 

Attn:  Code  0^-098  B 

Strategic  Systems  Project  Office 
Department  of  the  Navy 
Washington,  DC  20376 
Attti:  NSP-200 

Naval  Air  Systems  Command 
Department  of  the  Navy 
Washington,  DC  20361 

Attn:  Code  5302  (Aerospace  and  Structures) 
604  (Technical  Library) 

320B  (Structures) 


Naval  Air  Development  Center 
Director,  Aerospace  Mechanics 
Warminster,  PA  18974 

U.S.  Naval  Academy 
Engineering  Department 
Annapolis,  MD  21402 


Commanding  Officer  and  Director 
David  W.  Taylor  Naval  Ship 

Research  and  Development  Center 
Bethesda,  KD  20084 
Attn:  Code  042 
17 
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Naval  Underwater  Systems  Center 
Newport,  RI  02840 
Attn:  Dr.  R.  Trainor 


Naval  Surface  Weapons  Center 
Dahlgren  Laboratory 
Dahlgren,  VA  22448 
Attn : Code  G04 
G20 


Naval  Facilities  Engineering  Command  Attn:  Code  G04 

200  Stovall  Street  G20 

Alexandria,  VA  22332 

Attn:  Code  03  (Research  and  Development)  Technical  Director 

04B  Mare  Island  Naval  Shipyard 


14114  (Technical  Library) 

Naval  Sea  Systems  Command 
Department  of  the  Navy 
Washington,  DC  20362 

Attn:  Code  03  (Research  and  Technology) 
037  Ship  Silencing  Division) 
035  (Mechanics  and  Materials) 

Naval  Ship  Engineering  Center 
Department  of  the  Navy 
Washington,  DC  20362 
Attn:  Code  6105G 
6114 
6120D 
6128 
6129 


Vallejo,  CA  94592 

U.S.  Naval  Postgraduate  School 

Library 

Code  0384 

Monterey,  CA  93940 

Webb  Institute  of  Naval  Architecture 
Attn:  Librarian 
Crescent  Beach  Road,  Glen  Cove 
Long  Island,  NY  11542 
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Army 

Commanding  Officer  (2) 

U.S.  Army  Research  Office 
P.O.  Box  12211 

Research  Triangle  Park,  NC  27709 
Attn:  Mr.  J.  J.  Murray, 
CRD-AA-IP 

Watervllet  Arsenal 
MAGGS  Research  Center 
Watervllet,  NY  12189 

Attn:  Director  of  Research 

U.S.  Army  Materials  and  Mechanics 
Research  Center 
Watertown,  MA  02172 

Attn:  Dr.  R.  Shea,  DRXMR-T 

U.S.  Army  Missile  Research  and 
Development  Center 
Redstone  Scientific  Information 
Center 

Chief,  Document  Section 
Redstone  Arsenal,  AL  35809 

Army  Research  and  Development 
Center 

Fort  Bel voir,  VA  22060 


Air  Force 
Commander  WADD 

Wright-Patterson  Air  Force  Base 
Dayton,  OH  45433 
Attn:  Code  WWRMDD 

AFFDL  (FDDS) 

Structures  Division 
AFLC  (MCEEA) 

Cl.lef  Applied  Mechanics  Group 
U.S.  Air  Force  Institute  of  Technology 
Wright-Patterson  Air  Force  Base 
Dayton,  OH  45433 

Chief,  Civil  Engineering  Branch 
WLRC,  Research  Division 
Air  Force  Weapons  Laboratory 
Kirtland  Air  Force  Base 
Albuquerque,  NM  87117 

Air  Force  Office  of  Scientific  Research 
Bolling  Air  Force  Base 
Washington,  DC  20332 

Attn:  Mechanics  Division 

Department  of  the  Air  Force 
Air  University  Library 
Maxwell  Air  Force  Base 
Montgomery,  AL  36112 


NASA 

National  Aeronautics  and  Space  Administration 

Structures  Research  Division 

Langley  Research  Center 

Langley  Station 

Hampton,  VA  23365 


National  Aeronautics  and  Space  Administration 
Associate  Administrator  for  Advanced 
Research  and  Technology 
Washington,  DC  20546 

Scientific  and  Technical  Information  Facility 
NASA  Representative  (S-AK/DL) 

P.O.  Box  5700 
Bethesda,  MD  20014 
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other  Government  Activities 
Commandant 

Chief,  Testing  and  Development  Division 
U.S.  Coast  Guard 
1300  E Street,  NW 
Washington,  DC  20226 

Technical  Director 
Marine  Corps  Development 
and  Education  Command 
Quantico,  VA  22134 

Director  Defense  Research 
and  Engineering 
Technical  Library 
Room  3C128 
The  Pentagon 
Washington,  DC  20301 

Director 

National  Bureau  of  Standards 
Washington,  DC  20034 

Attn:  Mr.  B.  L.  Wilson,  EM  219 

Dr.  M.  Gaus 

National  Science  Foundation 
Environmental  Research  Division 
Washington,  DC  20550 

Library  of  Congress 

Science  and  Technology  Division 

Washington,  DC  20540 

Director 

Defense  Nuclear  Agency 
Washington,  DC  20305 
Attn:  SPSS 

Mr.  Jerome  Persh 
Staff  Specialist  for  Materials 
and  Structures 
OUSDR&E,  The  Pentagon 
Room  301089 
Washington,  DC  20301 

Chief,  Airframe  and  Equipment  Branch 
FS-120 

Office  of  Flight  Standards 
Federal  Aviation  Agency 
Washington,  DC  20553 


National  Acaden\y  of  Sciences 
National  Research  Council 
Ship  Hull  Research  Coimilttee 
2101  Constitution  Avenue 
Washington,  DC  20418 
Attn:  Mr.  A.  R.  Lytle 

National  Science  Foundation 
Engineering  Mechanics  Section 
Division  of  Engineering 
Washington,  DC  20550 

Picatinny  Arsenal 

Plastics  Technical  Evaluation  Center 
Attn:  Technical  Information  Section 
Dover,  NJ  07801 

Maritime  Administration 
Office  of  Maritime  Technology 
14th  and  Constitution  Ave.,  NW 
Washington,  DC  20230 

Maritime  Administration 
Office  of  Ship  Construction 
14th  and  Constitution  Ave.,  NW 
Washington,  DC  20230 
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Until  recently  only  the  theoretical  aspects  of  formulating  the  analysis 
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Proper  design  and  Implementatlolikof  user-oriented  computer  software  for 
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are  again  considered,  however,  primary  emphasis  has  been  placed  on  design 
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'^The  governing  equations  of  nonlinear  continuum  mechanics,  usually 
expressed  in  tensor  form,  are  given  in  matrix  notation  familiar  to 
most  structural  engineers.  Both  geometric  and  material  nonlinearities 
are  considered. 


The  classical  equations  of  the  Lagrangian  description  are  cast  into 
the  finite  element  method.  Specific  matrices  are  given  for  large 
deformations  of  a three  dimensional  solid  to  illustrate  the  procedure. 


A general  purpose,  user-oriented  software  system,  FINITE,  for  static 
linear  and  nonlinear  analysis  is  described.  FINITE  relies  upon  the 
POLO  supervisor  for  problem-oriented- language  translation,  data 
management,  and  dynamic  memory  allocation.  FINITE  supports  multi-level 
user-defined  substructuring  and  condensation  for  linear  and  nonlinear 
analysis.  Isolation  of  the  element  and  material  model  libraries  from 
the  system  enables  rapid  entry  of  new  modeling  capabilities., 


A number  of  nonlinear  example  analyses  are  presented  to  demons 
features  of  the  FINITE  system. 
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